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This thesis details the system identification and initial system validation of the an 
Ultra-Quiet Vibration Isolation Platform (UQP). With the move toward lighter and more 
flexible spacecraft, the effects of vibration are of immense concern. As natural or passive 
damping becomes less effective in controlling undesired vibrations, active vibration 
control becomes essential. The UQP uses a special configuration of the six degree of 
freedom Stewart Platform with piezoceramic strut actuators and geophone sensors. This 
combination gives an extremely sensitive and responsive six degree-of-freedom active 
vibration control system. Each actuator was designed to be controlled independently 
without coupling with other actuators. In order to develop control laws, the plant must be 
identified in terms of system zeros and poles and the uncoupled design validated. 
Dynamic modeling using parametric estimation methods can accurately describe a 
complex system. Using parameter estimation methods, models of the actuator system 
dynamics were obtained. A simple lead-lag controller was applied to individual actuators 
then all six actuators acting simultaneously to verify system coupling. Significant 
interaction between base adjoining actuators was discovered. 
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I. INTRODUCTION 



A. PURPOSE FOR RESEARCH 

With increased use of lightweight materials, spacecraft will become more flexible 
and susceptible to excitation from sources of vibration. With loss of stiffness and rigidity, 
passive damping alone may be insufficient to attenuate disturbances to an acceptable 
level. Active Vibration Control (AVC) is rapdily emerging as a method of eliminating or 
reducing unwanted vibration. The purpose of this thesis is to provide a detailed 
description and system identification of the Ultra Quiet Vibration Isolation Platform 
(UQP). It will be the basis for future research in the area of AVC. The UQP is an 
excellent sytem for AVC and can be applied to a broad range of disturbances. A 
disadvantage of the UQP is the complexity of the plant dynamics and kinematics. 
Accordingly the complexity of the control system design is increased. This thesis will 
charaterize the plant for future control design on the UQP. 

B. SCOPE 

The “plant” is the UQP which consists of six piezoceramic control actuators. This 
thesis will use parameter estimation methods to create a dynamic model of the UQP 
control actuators with the ultimate goal of extracting actuator zeros, poles and transfer 
functions. These actuators are designed to operate independently. Theoretically, each 
actuator can be controlled without coupling or interaction with the other actuators. This 
assumption is of extreme importance as it critically impacts upon the performance of the 
system and control system design. This assumption will be checked for validity by 
applying a lead-lag feedback compensator to an actuator and then to all six actuators 
simultaneously. 
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II. EXPERIMENTAL SETUP 



A. BACKGROUND 

In many spacebome applications the dynamics and control of vibrations must be 
addressed as a multiple degree-of-freedom (DOF) problem. Translations and rotations 
about all three axis must be considered. The Stewart Platform is the ideal mechanism for 
multiple DOF vibration control applications. 

1. The Stewart Platform 

The original motivation put forth by Stewart [Ref. 1] was to design a 
mechanism capable of simulating flight conditions to train pilots. In order to be realistic 
it had to be capable of translating and rotating in three directions just as a real aircraft. 
The original configuration consisted of a triangular plate and a rigid parallel base 
connected by six legs in a “linear coordinate” leg system. At each triangle vertex two 
legs were attached in a three DOF joint mechanism. These legs, with hydraulic actuators, 
were mounted to two DOF joints at the base. The advantages of choosing this 
configuration were inherent rigidity and absence of bending moments. Additionally with 
this configuration the only forces present were in the plane of the leg. A similiar six leg 
arrangement had also been used by a machine designed to study tire to ground forces [Ref 
2], In this system computer controlled jacks were used as actuators. A cubic 
arrangement was devised such that the relationship of one actuator to any other is the 
same (see Figure 2.1). This arrangement showed inherent stability and the capability for 
accepting large moments. 
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Figure 2.1 
Cubic Arrangement 

Stewart Platform mechanisms have been the subject of studies as multiple DOF parallel 
systems [Ref. 3 and Ref. 4:p 46]. The Stewart Platform configuration has shown wide 
applicability from motion simulators to robotics and now Active Vibration Control. The 
UQP employs the cubic configuration of the Stewart Platform which has an extremely 
important property with respect to AVC. With the exception of inertial loads and gravity 
forces, all other forces are carried axially. The significance is that if axial forces due to 
vibration can be eliminated the vibration is eliminated [Ref. 4:p. 46]. 

B. HARDWARE CONFIGURATION 



The experimental hardware is divided into two sections: the UQP and the Power, 
Control and Analysis Hardware. 
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1 . 



Ultra-Quiet Vibration Isolation Platform (UQP) 



The UQP is an adaptation of the Stewart Platform for AVC designed and 
built by CSA Engineering INC. It consists of a circular plate attached to a rigid base by 
six variable length actuators (see Figure 2.2). 




Figure 2.2 

Ultra-Quiet Vibration Isolation Platform 

The base consists of an aluminum plate supported by aluminum stringers and longerons. 
The base is used to simulate a rigid spacecraft. Mounted inside one of the “spacecraft” 
bays is an AURA Bass Shaker which is used as a known disturbance source. The entire 
aparatus is bolted to a 38001b NEWPORT optical table. 

a) Control Actuator Struts 

The UQP uses six mutually orthoganol struts to provide control 
over six degrees of freedom. Each strut consists of a piezoceramic actuator and a 
geophone sensor. The use of piezoceramics for shape and vibration control is becoming 
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increasingly widespread. For interested readers, further information on piezoceramics 
and their applications in vibration control can be found in Reference 8. By taking 
advantage of the cubic configuration, all six struts can be considered linear motion 
actuators [Ref 1]. Stewart Platform based AVC mechanisms using the same cubic 
configuration with Terfanol-D actuators [Ref. 4] and voice coil actuators [Ref. 5] have 
also been developed. These platforms displayed significant reduction of vibration over 
passive means alone. Figure 2.3 is a basic diagram of each actuator. The actuator stroke 
is approximately 50 microns. In the passive (open loop) mode it provides moderate 
damping to low frequency vibrations. In the active (closed loop) mode it provides 
damping to higher frequency vibrations. Passive damping is provided by a flexure with 
damping material. An extremly sensitive geophone measures velocity. This velocity is 
the error signal fed into the control law. 




Figure 2.3 
Actuator 
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In general applications, geophones consist of wire coils supported by soft springs under 
the influence of a magnetic field. Vibrations cause movement of the magnet relative to 
the stationary coil inducing a voltage proportional to velocity [Ref. 6:p 449]. The 
stiffness K refers to the stiffness of the passive stage. The stiffness K a corresponds to the 
piezoceramic stack actuator (PZT). 



b) Disturbance Shaker 

To create a disturbance source against which the performance of 
the UQP and applied control can be measured, a known disturbance source was mounted 
to the base. A sinusoidal waveform from a signal generator is amplified and input to a 
AURA Bass Shaker (see Figure 2.4). This simulates a cyclic disturbance. It has a 
resonant frequency at 42 Hz which is used as the disturbance frequency in initial 
experiments. 



Disturbance 

Shaker 




Figure 2.4 
Disturbance Shaker 
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2 . 



Power, Control and Analysis 



To operate the UQP and extract meaningful experimental data requires 
several important components (see Figures 2.5 & 2.6). 



Control Input-Sensor Output 



KEPCO Operational 
Amp. & Pwr. Supply 




PWR Supply 






AVCS CX- 



ft It 



HP 33120A Function/ 

Arb. Waveform Generator 



Sensor Output- Channel of Interest 








PC 


dSPACE 







Figure 2.5. 

Power, Control and Analysis Hardware Configuration 




Figure 2.6 

High Voltage Power Supply (HVPS)(top) and Active Vibration Control System (AVCS) 
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The current configuration requires a MATLAB capable PC with an RS232 I/O port. The 
control design is done in MATLAB. Once the design is completed, the MATLAB code is 
translated into C code. The C code is converted to machine language and down loaded 
via the I/O port to the CSA Active Vibration Control System (AVCS) (see Figure 2.6). It 
provides the interface for control implementation and geophone sensor output. Its main 
component is a digital signal processor (DSP). Once the machine code has been loaded 
the AVCS applies the control to the actuators. It receives the feedback signal from the 
geophones, processes and feeds it into the control algorithm. The output is sent to a HP 
35665A Dynamic Signal Analyzer for analysis. The CSA High Voltage Power Supply 
(Figure 2.6) supports the power equirements of the piezoceramic actuators and is not 
required for use in the open loop mode. 
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III. SYSTEM IDENTIFICATION 



In order to design an effective compensator or controller for a dynamic system the 
plant must be well known or characterized. By the use of system identification 
techniques a plant model can be obtained. The use of a model is extremely beneficial in 
cases when the dynamics of a plant are not well known or complex. Although articles 
have been written investigating the dynamic and kinematic complexities of Stewart type 
platforms [Ref. 3.], active vibration control using Stewart Platform based mechanisms is 
a relatively new field. Experimental dynamic modeling can be used to bypass the need 
for a complete theoretic dynamic analysis of the UQP. The goal of identifying the UQP 
plant was to obtain an accurate plant transfer function which can assist control design. Of 
particular interest were system modes or resonances and interaction or couplings among 
actuators. An initial key assumption on the actuators is that they each act independently 
(uncoupled) and can be controlled as such. Validation of this assumption is one of the 
goals of identification of each plant 

A. DYNAMIC MODELING 

There are two categories of system identification used in dynamic modeling: Non- 
parametric and Parametric. Non-parametric methods are used to estimate transient 
response, spectra and frequency functions in order to gain basic knowledge of the system. 
Parametric methods are used to obtain the mathematical parameters or elements used by 
well know system modeling algorithms. These models are broken down into two 
subcatagories. “Tailor Made” models are built based on physical principles (i.e. the 
parameters have some type of physical interpretation). “Ready Made” or “Black Box” 
models are general in nature and use parameters that have no direct physical 
interpretation. They are used to characterize the relationship of output to input of a 
system. By treating the UQP as a “Black Box” and analyzing the input-output 
relationship it is possible to obtain an accurate model of the system by using a Ready 
Made model. [Ref. 10] 
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1. Background 



The concept behind the use of transfer functions is to describe the 
relationship between the input and output of a system and can be expressed in terms of 
either continuous or discrete time. The use of discrete time representations are preferable 
when the data is sampled in nature. 



u(t) 





Transfer 






Function 





-►y(t) 



Figure 3.1 

Input-Output Relationship 



One of the most important tools used to gain an understanding of the behavior of a system 
is by looking at its impulse response. For a linear system (plant) the impulse response is 
the output when the applied input is an impulse at t=0. 



5(0 






Plant 



► y(t)=h(t) 



Figure 3.2 

Impulse Response Relationship 



If the impulse response h(t) of a system is known, the output to a given input can readily 
be determined by convolution of the input with the impulse response: 

y(t) = u(t)*h(t) (3.1) 



or: 

y(t) = Y.h{t)*u{t-m) (3.2) 

The impulse response is a key element of non-parametric system identification. 

An essential element of system identification is characterizing the output. 
In an ideal system the output would simply be the input modified by the system transfer 
function. In physical systems however, the output also contains elements resulting from 
internal and external disturbances (see Figure 3.3). 
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Disturbance 



Disturbance 



U(t) 



4 > 



>- 



Plant=h(t) 




y(t)=y(t)+G(t) 



Where: 

y(t)=Total Output 
\\f(t)=Undisturbed Output 
a(t ^Disturbance term 



Figure 3.3 

In order to develop a model of the input output relationship, define: 

y(t) = y/(t) + a(t) 
y/(t) = G(q,p)u(t) 

o(t) = H(q,p)e(t) 
where £(t) is white noise, p is the parameter vector: 

( Pl> 

p= i 



\PnJ 



and q T is the shift operator based on sampling interval T: 

q T y(0=y(t+T) 

<1t y( t ) = y(t-T) 



(3.3 a,b,c) 



(3.4) 



(3.5 a,b) 



Introduce the general difference equation: 

y(t) + a^y{t — \T) + a 2 y{t -2T)+. ..+a N y(t - NT) = b 0 u(t) + b x u(t - lT)+...+b L u(t - LT) 

(3.6) 



where T is the sampling interval of the discrete time representation. This standard 
expression is a specialized geometric series that relates input and output of a system using 
past values of input, output and a set of parameters a, and bj. By applying the shift 
operator q T> equation (3.6) becomes: 
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(3.7) 



y(t)[l+aiq T ] + a 2qT 2 +'"+ a n4T n '\ = u ( t )\bo+b\qx l +b2qT 2 +---+b n q T n ^ 

y(t ) _ [^0 - > rb l q T l J rb 1 qj- +• • -+b n q T n ] 
Jl+a 1 < 7 7 - 1 +a 2 <? 7 - 2 +"-+^ 7 ’ n ] 

Define: 

G(-7»P) = 2 7T (3-8) 

u(t) 



and introduce a time delay of nk samples on the input: 

B(q) 



G(q,p) = 



- nk +b l q ( 7 r nk - l) +b 2 q ( f nk - 2) +- 



■+b n q<f nk - nb) ] 



A{q) 



[l+^l^r 1 +a 2 qT 2 +' • ' +a na^T na ] 



(3-9) 



Changing variables to be consistent with the notation used by Ljung and Glad [ref. 11] 
and the MATLAB System Identification Toolbox (SITB) [Ref. 12]: 



[ 1+ /l < 77' 1+ /2 < 7r 2+ ‘" + /n/^7’^] 



(3.10) 



Applying to equation (3.3b) and equation (3.6): 
yf(t) + / 1 y/(t - T)+- • •+ f n j y/(t - njT ) = bQu(t- nkT )+ • • •+ b n b u(t-(nb + nk)T (3.11) 



Similarly the disturbance term becomes: 



H(q,p) = 



C(q) _ [* +c l Vi 1 + c 2 <lT 2 +• - ‘ +c nc<lT nC ] 

D(q) ^+d\qj l +d2qT 2+ "-+d n d^T nCl \ 



(3.12) 



The coefficients bj, Cj, d, and f comprise the parameter vector p. These coefficients 
represent the unknown system parameters. Essentially they approximate a complex 
system’s physical parameters but have no direct correlation to any specific physical 
quantity. The parameters nb, nc, nd and nf characterize the order and type of ready made 
model. The model constructed by equations (3.3) through (3.12) is shown in equation 
(3.13). It represents the Box-Jenkins model 1 [Ref. 10]. 

y(t) = G(q,p)u(t) + H(q,p)e(t) (3.13) 



1 Named after the statisticians G.E.P. Box and G. M. Jenkins 
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The following models are variations of the Box-Jenkins model. The most 
simple one is the case where the disturbance is not modeled. 

y(t) = G(q,p)u(t) + e(t) (3.14) 

Equation (3.14) is known as the Output Error method. The difference between the actual 
output and undisturbed output is manifested in the white noise term £(t). The next model 
utilizes the same poles for the input and disturbance (noise) models defined as: 

F{q) = D{q) = A(q) = 1 + a, q~ T ' +--a w qT (3-15) 



Applied to equations (3.10), (3.12) and (3.13) yields: 



W) A{q) UW+ A{q) 



£(0 



(3.16) 



or: 

A(q)y(t) = B(q)u(t)+C(q)e(t) (3.17) 

This is known as the ARMAX (AutoRegression Moving Average eXtra input) model. 
The final model that will be described is the AutoRegression eXtra input(ARX) model 
(equation 3.18). [Ref. 10] 

A(q)y(t)=B(q)u(t) + £(t) (3.18) 

The next step is to describe how these ready made techniques model a 
system based on past values of output and input. Essentially these algorithms attempt to 
predict the output based on a given set of input-output data. Introduce the concept of 
One-Step- Ahead prediction of output. From equations (3.6) and (3.14), we have: 
y(t) = -a x y{t- X)-..-a im y(t — na) + b x u(t - nty+.^+b^uit -nk — nb + 1) + £(t) (3.19) 



Since e(t) is white noise and cannot be predicted, the prediction becomes: 

y(t, p ) = -a x y(t - 1)-. . -a w y(t - na) + b, u(t - nk)+. . .+b nh u(t - nk-nb + l ) (3.20) 

Expanding to the general case of equation (3.13) and dividing out the noise transfer 



function: 

H~' (q, p)y(t) = H~' (q, p)G(q,p)u(t ) + £(t) (3.21 ) 

y(t) - y(t) + //-' (q, p)y(t) = H' x (<?, p)G(q, p)u(t) + £(t) (3.22) 

y(t) + [- 1 + H~' (q, p)]y(t) = H~' (q, p)G(q, p)u(t) + £{t) (3.23) 

y(t) = [\-H-\q, p)]y(t) + H-\q,p)G{q,p)u{t) + £(t) (3 .24) 
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the prediction becomes: 

y(t, p) = [1 - H~ l (q, p)]y(0 + H~\q, p)G(q, p)u{t) (3.25) 

Since equation (3.25) only contains past values of the output y(t) and input u(t) the 
difference between the prediction and the output (the prediction error) is: 

E(t,p) = y(t)-y(t,p) = e(t) (3.26) 

where £(t) is white noise. [Ref. 9: p. 56 and Ref. 10: p.235] 



B. FREQUENCY RESPONSE OF ACTUATORS 



The diagram in Figure 3.4 is a modification of the experimental set up in Figure 
2.5 used to extract the input-output data. 
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Figure 3.4 
Experimental Setup 



The purpose behind measuring the frequency response of each actuator is to obtain both 
basic insight into the plant and a model verification baseline. Each actuator was excited 
by a 50 mV “pink noise” source provided by the HP Dynamic Signal Analyzer (HPDSA). 
The “pink noise” is random noise which provides 3dB roll off per octave. The motivation 
to use pink noise is to place an equal amount of energy in each octave band. The 
response is sensed by the geophone of the excited channel (actuator) and fed into the 
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HPDSA. The HPDSA computes and plots the time averaged frequency response which 
are shown in Figures 3.5 through 3.10. 



Strut #1 Magnitude 





Figure 3.5 

Frequency Response of Actuator #1 



Strut #2 Magnitude 





Figure 3.6 

Frequency Response of Actuator #2 
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Magnitude(dB) . , Deg Magnitude(dB) 



Strut #3 Magnitude 




Figure 3.7 

Frequency Response of Actuator #3 

Strut #4 Magnitude 





Figure 3.8 

Frequency Response of Actuator #4 
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Start #5 Magnitude 





Figure 3.9 

Frequency Response of Actuator #5 



Strut #6 Magnitude 





Figure 3.10 

Frequency Response of Actuator #6 
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The plots show that the plant transfer functions for the actuators are basically the same. 
The resonant dynamics of the struts are readily discernible, specifically the actuator 
natural frequency at 1 .4 kHz. 



C. MODEL STRUCTURE SELECTION 

The type of system and method of data collection are major factors determining a 
suitable model type for a given application. For example, the Output Error, OE, would be 
best applied when the disturbances (which are not modeled) are small. The ARX is the 
most widely used because it is the simplest and serves as a good baseline to select a 
model. However, it uses the same poles to model noise as those used to model the 
system. This can create problems in system identification requiring the use of higher 
order models than would otherwise be necessary. The Box-Jenkins method models the 
system dynamics and disturbances separately with no common parameters. The ARMAX 
also models the disturbance but uses the same poles as the system dynamic model. It 
differs from the ARX by the addition of the C(q) polynomial. The Box-Jenkins and 
ARMAX both provide detailed system models and can be extremely accurate. The 
decision on which to use is based on the nature of the disturbance. If “noise” enters the 
process in the early stages or is carried through the plant (see Figure 3.3) the ARMAX is 
preferable over the Box-Jenkins model which is more adept to characterizing noises that 
enter later in the system. In the case of the UQP, the ARMAX would appear to be the 
preferred method if there is interaction between actuators because the resulting 
disturbances would enter into the process early on. Additionally the disturbance caused by 
another actuator should show similar dynamics as the plant under experimentation. The 
flow chart in Figure 3. 1 1 represents the process used in building a model of each actuator. 
As Figure 3.1 1 shows, the process of system identification and modeling is a step-by-step 
iterative process. The basic idea was to develop and validate a basic ARX model structure 
(i.e. na, nb and nk) for all six actuators and then apply it to the more complex ARMAX 
and Box-Jenkins methods and determine the best model. The application of this process 
on actuator number one will be detailed in the following sections. [Ref. 10] 
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Figure 3.1 1 

Model Selection Flow Chart 

The computer program used to accomplish this can be found in appendix A. 

1. Data Collection 

The process of system identification and dynamic modeling relies on a set 
of plant input-output data. For the UQP each actuator was provided an excitation and the 
response of all six actuators (as measured by each respective geophone) recorded (see 
Figure 3.4). Using the HPDSA as the excitation source, 50mVRMS “Pink Noise” was 
applied to the actuator via the 25 pin connection on the front of the AVCS (see Figures 
2.6 and 3.3). Under normal operation the ribbon cable on the front of the AVCS is 
connected to close the control loop. By disconnecting the cable, each actuator can be 
accessed individually. The geophone sensor output was connected via a BNC interface 
box to a dSPACE system. The dSPACE was used to gather, display and save all six 
channels of data simultaneously. For comparison purposes the active channel sensor 
output (actuator under excitation) and excitation source was also input to the HPDSA. 
For each actuator the output corresponding to the excitation source was measured at a 
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sampling frequency of 10 kHz for two seconds. This gave 20,000 input/output sample 
points. The choice of sampling frequency was based on the desire to obtain a sufficient 
number of sample points from which to construct and validate a model. The first 10,000 
(input and output) sample points are used by the parameter estimation algorithms to build 
the model. The last 10,000 points (of input) are used in model validation. If the sample 
points used to build the model are also used in validation, the true accuracy and validity 
of the model will be corrupted. 



a) Data Preparation 

Prior to using the data to build a model it is necessary to pre-treat 
the data. The most common factors adversely affecting collected data are: 

1 . Drift, offset, trends and low frequency and/or periodic disturbances. 

2. Outliers or faulty data points. 

3. High frequency disturbances. 

Before treating for these possible problems the data is given a preliminary analysis to see 
if they do in fact exist or could be the possible source of erroneous results [Ref. 9:p 386]. 

External sources are the main cause of drifts, trends and low 
frequency disturbances that are either impossible or undesirable to model. This can be 
remedied by removing external disturbances, using the noise model to account for the 
disturbances and/or using an algorithm to de-trend the data. Due to the extreme 
sensitivity of the geophone sensors it is impossible to remove all the external 
disturbances. The latter two options were selected. It is assumed that the external 
disturbances received by the geophones are relatively small compared to those entering 
into the process early on and will be adequately addressed by the noise model. The offset 
problem is a result of the failure of input-output data to correspond in a consistent 
manner. This causes the model to waste parameters in attempt to adjust these levels. 
This is compensated by a MATLAB function included in the computer code (see 
appendix A.). [Ref. 12:p 1-63] 
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In any data collection effort there are obviously erroneous values. 
However, removal of these “outliers” requires caution. Removal can either be manual 
(i.e. selecting bad values by hand) or by some type of automatic algorithm. The best way 
to determine if outliers are problematic is to check the model residuals for excessively 
large values that stand out. This method was chosen due to the vast amounts of sampled 
data used. [Ref. 12:pl-63] 

The final source of data deficiency is the presence of high 
frequency disturbances outside the region of interest. In this case the primary concern is 
the performance of the actuators up to and including the resonance frequency. A 10th 
order Butterworth filter with a pass band below approximately 1.5kHz was used to 
eliminate disturbances outside the range of interest [Ref. 1 1: p 272]. 

2. Delay and Parameter Number/Structure Selection 

Prior to applying one of the aforementioned parametric models to the 
gathered data, the input delay nk, and the structural elements na, nb, nc, nd and nf must 
be determined. This was done using the Model Structure Selection functions in the 
MATLAB SITB. The SITB has functions to calculate ARX-based and OE-based 
structures. Figure 3.12 represents a flow chart of the algorithm used to select the delay 
and structure. Since the assumption was that the ARMAX model would provide the best 
model for the UQP, an ARX-based structure selection algorithm was applied (see 
appendix A). 
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Figure 3.12 

Structure Selection Algorithm 

The process of selecting na and nb as well as the number of total parameters from which 
they are derived is an iterative process. Since the input delay should be the same no 
matter what the order of the model is, a second order ARX model will serve as the 
foundation of this process. Table 3.1 shows the results of nk values ranging from 1 to 10 
applied to a second order ARX model. 
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Act. 


nk=l 


nk=2 


nk=3 


nk=4 


nk=5 


nk=6 


nk=7 


nk=8 


nk=9 


nk=10 


i 


-5.5975 


-5.6122 


-5.4704 


-5.4708 


-5.4659 


-5.4686 


-5.4755 


-5.5046 


-5.5184 


-5.5077 


2 


-5.3380 


-5.3175 


-5.2293 


-5.2301 


-5.2297 


-5.2375 


-5.2463 


-5.2602 


-5.2596 


-5.2426 


3 


-4.6214 


-4.4221 


-4.4132 


-4.4391 


-4.4432 


-4.4293 


-4.4357 


-4.4282 


-4.4252 


-4.4256 


4 


-5.5237 


-5.5711 


-5.3791 


-5.3872 


-5.3928 


-5.3946 


-5.4007 


-5.4106 


-5.4054 


-5.3905 


5 


-5.0731 


-5.1164 


-5.0082 


-5.0077 


-5.0044 


-5.0028 


-5.0089 


-5.0239 


-5.0259 


-5.0130 


6 


-5.8464 


-5.8072 


-5.7362 


-5.7301 


-5.7204 


-5.7390 


-5.7443 


-5.7598 


-5.7747 


-5.7610 



Table 3.1 
Input Delays 

The numbers represent the logarithms of a quadratic loss function based on the selected 
nk for na=nb = 2 (second order model). Once an input delay has been selected a plot of 
the number of parameters used versus loss function is produced. Based on the results in 
Table 3.1, the best choice for an input delay for actuator number one appears to be nk= 2. 
Figure 3.13. is the resulting plot. 
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Figure 3.13 

Parameters Vs. Loss Function for Actuator #1 nk = 2 

The multiple results for a given number of parameters is due to the fact that there can be 
several different structures (i.e. combinations of na and nb ). The SITB has functions that 
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automatically compute the structures based on the minimization of Akaike’s Information 
Theoretic Criterion (AIC) and Rissanen’s Minimum Description Length (MDL)[Ref. 1 1], 
After a structure has been selected it can be applied to a parametric estimation method to 
develop a model. Even with the availability of functions capable of selecting the 
optimum number of parameters based on minimization techniques, it can be seen in 
Figure 3.13 that there is a point of diminishing returns beyond which the addition of more 
parameters is of little benefit. Further, the addition of more parameters than necessary 
will actually cause a loss in model accuracy. Referring to Figure 3.13 the structures 
chosen are listed in table 3.2. 



Number of Parameters 


Structure [na nb nk] 


Notes 


30 


15 15 2 


AIC Based 


22 


14 8 2 


MDL Based 


18 


14 4 2 


Automatic 



Table 3.2 

Structures For Actuator #1 Input Delay nk= 2 

The notes column provides information on how the structure was obtained. AIC and 
MDL based structures were automatically generated along with the plot in Figure 3.13. 
The term “Automatic” means the number of parameters in the first column were manually 
extracted from the plot in Figure 3.13 and entered in to an SITB function that 
“Automatically” selects the structure. The term “Manual” refers to the complete structure 
being manually selected. 

3. Model Structure Application and Evaluation 

After determining the ARX model structure for actuator number one, a 
model was computed using the algorithm in appendix A. There are numerous methods 
for measuring the accuracy and validity of a model. The ones utilized to validate the 
actuator models are described as they were applied to the ARX model of actuator number 
one using the structure [na=15 nb=15 nk=2]. 
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models are described as they were applied to the ARX model of actuator number one 
using the structure [na=15 nb=15 nk=2], 

a) Frequency Response 

One of the quickest ways to determine a model’s validity is to 
compare the Bode plot of the system to that of the model. 




Frequency (Hz) 



(a) Magnitude 




Frequency (H2) 

(b) Phase 
Figure 3.14 

Frequency Response of ARX[15 15 2] Model (Actuator #1) 
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The model’s frequency response is similar to that of the actual system (see Figure 3.14). 
The areas of critical concern are the model behavior at resonance and zero dB crossing for 
magnitude and -180 degree crossing for phase. The phase is reasonably close however, 
the magnitude is significantly different. The general shape and resonance peak are in 
agreement but there is a distinct 10-20dB difference between the actuator and the model. 

b) Output Comparison 

In evaluating a model, as stated earlier, it is important to use 
different input sample data to insure an accurate assessment of the model. In modeling the 
actuators a 10,000 sample model validation set was held aside for this puppies. Applying 
the validation input to both the model and the system and comparing the outputs, gives an 
accurate measure of the model’s validity and how close it truly simulates the actual 
system. Figure 3.15 compares the model output to the actuator output in response to the 
same input. 



Plot of ARX MODEL [15 15 2] Output vs. Actual output 




Sample 



Figure 3.15 

Actual Vs. Model Output For ARX [15 15 2] (Actuator #1) 
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c) Auto and Cross Correlation Functions 



Figure 3.16 provides plots of the auto correlation function of 
residuals and cross-correlation function between the input and residuals from the output. 
The difference between the model output and actual output is called the residual. 
Basically, it is what is left unaccounted for by the model. Recalling equation (3.26): 

E(t, p ) = y(t) - y(t, p) = £(t) 

where e(t) represents the residuals. If £(t) is purely “white noise” it is independent of the 
input. If there is correlation between £(t) and u(t), there are elements in the output from 
the input that are not explained by the model. 



Correlation function of residuals. Output# 1 




- 0 . 5 1 1 1 1 1 1 

0 5 10 15 20 25 

Cross corr. function between input 1 and residuals from output 1 




Figure 3.16 

Auto and Cross-Correlation Functions for ARX [15 15 2] Model for Actuator #1 

The horizontal lines in Figure 3.16 indicate a 99% confidence level. If the function 
crosses these lines, there is a correlation between s(t) and u(t) at that point. From Figure 
3.16 there appears to be correlation between £(t) and u(t-l). This indicates that the input 
delay might need to be modified to one. The auto-correlation function is acceptable. A 
plot of residuals versus sample number for 100 samples is given in Figure 3.17. 
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0.1 



Residual of ARX Model [15 15 2] 
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Sample 

Figure 3.17 

Residuals Vs. Sample ARX[15 15 2] For Actuator #1 

Over this range of samples the residuals are relatively low. For an average output 
magnitude of approximately 0.18, the average residual magnitude is on the order of 0.03 
or roughly 17% of the output magnitude. Although this percentage is a very rough way of 
comparison and cannot be considered by itself to determine model validity, it does give a 
basic idea of how much is left unexplained by the model. A much better method is to use 
the mean square fit. [Ref. 10] 



d) Zero-Pole Plot 

The final evaluation tool is the zero-pole plot (see Figure 3.18). 



30 



Zero Pole Plot ARX Model (1515 2] Model 




Figure 3 18 

Zero-Pole Plot of ARX [15 15 2] For Actuator #1 

The ellipses correspond to confidence limits to three standard deviations [Ref. 12], This 
plot shows that there are numerous possible zero-pole cancellations indicating that a 
model of lesser degree is desirable. 

4. Model Application and Evaluation (Data Prefdtered) 

Based on analysis the ARX [15 15 2] is not a very good model. Following 
the flow chart the same analysis will be repeated on data pretreated with a 10 th order 
Butterworth filter. 



a) Frequency Response 

By prefiltering the data the correlation between the model 
frequency response and the actual output frequency response is improved. However, 
there is still an approximate 1 0 dB difference between the two responses. 
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Blue=Actuator#1, Red=ARX M0DEL[15 15 2) *Data Prefiltered* 




(a) Magnitude 



B!ue=Actuator#1 , Red=ARX MODEL[15 15 2] ‘Data Prefiltered* 




Frequency (Hz) 



(b) Phase 
Figure 3.19 

Frequency Response of ARX[15 15 2] Using Prefiltered Data (Actuator#l) 

b) Output Comparison 

Figure 3.20 shows an improvement in the correlation between the 
model and actual output. 
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Plot of ARX MODEL [15 15 2] Output vs. Actual output ‘Data Prefiltered* 




Sample 



Figure 3.20 

Prefiltered ARX[15 15 2] Model Output Vs. Actual Output 



c) 



Auto Correlation and Cross Correlation Functions 



As a result of prefiltering there is a degradation of the auto and 
cross correlation functions (see Figure 3.21). However, the auto-correlation of residuals 
is relatively small For cross correlation functions where the crossing occurs in negative 
lag, this is an indication that the data was possibly collected during feed back vice a 
problem with the model [Ref. 10:p 284] 

Correlation function of residuals. Output # 1 




0 5 10 15 20 25 

Cross corr. function between input 1 and residuals from output 1 




Figure 3.21 

Auto and Cross Correlation Functions of Prefiltered ARX[15 15 2] (Actuator #1) 
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d) Residuals 



Comparing Figure 3.22 with Figure 3.20, the former residual plot 
is misleading. This anomaly occurs when large order models (large number of 
parameters) using prefiltered data are applied to the SFTB function. This function uses 
statistical means to calculate the prediction error and the residuals. A simple plot of the 
actual output minus the model output is shown in Figure 3.23. 




Sample 

Figure 3.22 

Residuals of Prefiltered ARX[15 15 2] (Actuator #1) 



Actual Output - Model Output ’Data Filtered* 




Sample 



Figure 3.23 

Actual Output Minus ARX[15 15 2] Output (Actuator#!) 
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The average “residual” magnitude for untreated and pretreated data are 0.0041 and 
0.0032 respectively, or 2% and 17% of the average output magnitude. As previously 
stated, these are very rough measures of the accuracy but there is a clear improvement in 
the model by prefiltering the data. 

e) Zero-Pole Plots 

The zero-pole plot of the prefiltered model also shows 
improvement over the unfiltered model (see Figures 3. 18 and 3.24). 



Zero Pole Plot ARX Model [15 15 2] Model 'Data Preliltered* 




Figure 3.24 

Zero-Pole Plot For Prefiltered ARX [15 15 2] (Actuator #1) 

However, there are still numerous possible zero-pole cancellations, indicating that use of a 
lower order model is necessary. 

5. Model Structure Acceptance for Actuator Number One 

The goal of the model structure application and evaluation process is to use 
the ARX based technique to find an acceptable model structure that can be applied to 
more accurate ARMAX and Box-Jenkins methods. Knowing that prefiltering data 
provides a better model, the process is repeated. Through the first iteration it was also 
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determined that a delay change would possibly be beneficial. Returning to the delay 
selection process, Figure 3.25 plots the loss function for a given number of parameters. 
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Figure 3.25 

Number of Parameters Vs. Loss Function For nk= 1 (Actuator #1) 
The structures evaluated are contained in table 3.3. 



Number of Parameters 


Structure [ na nb nk] 


Notes 


30 


15 15 1 


AIC/MDL Based 


24 


14 10 1 


Automatic 


20 


7 13 1 


Automatic 


18 


7 11 1 


Automatic 


15 


69 1 


Automatic 


13 


67 1 


Automatic 


13 


76 1 


Manual 



Table 3.3 

Structures For Actuator #1 Input Delay nk= 1 
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After a set of candidate structures was obtained they were used to generate ARX models. 
The process of model selection and evaluation outlined in Figure 3.11 and described in 
sections 3 and 4 was applied to each candidate structure based ARX model. The most 
suitable candidate found was ARX [na=7 nb=6 nk= 1] (see Figures 3.26 though 3.32). 
Table 3.3 lists the [na=7 nb = 6 n£=l](or [7 6 1]) structure as a manual selection. In the 
first iteration the [7 13 1] structure was selected. Upon repeating the selection process for 
actuator number two, the [7 6 1] structure was generated automatically by the algorithm 
(see appendix A) based on the choice of 13 parameters. The performance of the resulting 
model was judged superior and the [7 6 1] structure was applied to actuator number one 
where it also yielded the best ARX model (see Figure 3.32). However, two design trade- 
offs were made. First, the auto and cross correlation functions (see Figure 3.28) showed 
degraded performance, indicating there are still unaccounted input elements in the output. 
However, it does not affect the frequency response and other criteria used to validate the 
structure. The second trade-off was the zero-pole plot (see Figure 3.31 ). The close 
proximity of two poles and two zeros indicates that a lower order model would be 
desirable. However, this proved to not be the case. The use of a lower order model 
showed declining performance over the ARX model using the [7 6 1] structure. 
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Frequency (Hz) 
(a) Magnitude 




Frequency (Hz) 

(b) Phase 
Figure 3.26 

Frequency Response of Prefiltered ARX[7 6 1] (Actuator #1) 



38 



Plot of ARX MODEL [7 6 1] Output vs. Actual output 'Data Prefiltered* 




Sample 



Figure 3.27 

Actuator #1 Output Vs. ARX[7 6 1] Model 



Correlation function of residuals. Output# 1 




Cross co it. function between input 1 and residuals from output 1 




Figure 3.28 

Auto and Cross Correlation Functions of ARX[7 6 1] Model (Actuator #1) 
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Residual 



Residual of ARX Model [7 6 1] *Data Prefiltered* 




Sample 
Figure 3.29 

Residua] Vs. Sample for ARX[7 6 1] (Actuator #1) 



Actual Output-ARX[7 6 1] Model Output *Data Prefiltered* 




Sample 
Figure 3.30 

Actual Output and ARX[7 6 1] Model Output Difference (Actuator #1) 
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Zero Pole Plot ARX Model [7 6 1] Model ‘Data Prefiltered 
1 .5 , . . . ■ ■ 1 




Figure 3.31 

Zero-Pole Plot for ARX[7 6 1] Model (Actuator #1) 
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Strut #1 Magnitude 




10 3 



Frequency (Hz) 
(a) Magnitude 




Frequency (HZ) 

(b) Phase 

Figure 3.32 

ARX Model Frequency Response Comparison 
B lue= Actuator# 1, Red=[7 6 1], Magenta=[15 15 1], Black=[7 13 1] 
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By increasing the delay the correlation functions improve but performance in all other 
respects is degraded. The 99% confidence limit may also be too constrictive. Confidence 
limits of 95% have been used [ref. 10 p.447] to constrain correlation functions. These 
factors governed the decision to accept the degraded correlation functions as a design 
trade off. In comparison with models based on various structures (Figure 3.32) this 
model is acceptable. Table 3.4 summarizes numerical results. 



Measure 


Result 


Average Output Magnitude 


0.1870 


Average Residual Magnitude 


0.0064 


Percentage of Average Output 


3.4% 


Average Output Difference Magnitude 


0.0024 


Mean Square Fit 


0.1562 



Table 3.4 

Prefiltered ARX[7 6 1] Numerical Results (Actuator #1) 

This process was repeated to find the structures for an ARX based model for actuator 
numbers two through six. This process involved the analysis of hundreds of 
validation/verification plots and numerical results. Only the validation plots of the 
accepted structure for each respective actuator were included in this thesis and are shown 
in appendix B. 



6. Model Structure Acceptance for Actuator Number Two 

The process outlined in section five was repeated for actuator number two. 
The plots resulting from the acceptance process are listed in appendix B (Figures A.l 
through A.l). Based on a delay of nk=l from table 3.1 and resulting parameter number 
plot (see appendix B Figure A.l) the structures selected for evaluation are listed in table 
3.5. 
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Number of Parameters 


Structure [ na nb nk ] 


Notes 


26 


12 14 1 


AIC Based 


25 


11 14 1 


MDL Based 


24 


11 13 1 


Automatic 


22 


11 11 1 


Automatic 


18 


99 1 


Automatic 


15 


96 1 


Automatic 


13 


76 1 


Automatic 


12 


66 1 


Automatic 



Table 3.5 

Evaluated Model Structures for Actuator #2 

After evaluating the graphical and numerical results of the structures listed in table 3.5, 
the [7 6 1] structure was identified as the best structure. It provides excellent 
performance for a minimal number of zeros and poles. Figure 3.33 compares the [7 6 1] 
based model with the next best models. 
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Strut #2 Magnitude 
55 n 1 1 — i — i — r-rn 




Frequency (Hz) 



(a) Magnitude 




10 3 



Frequency (Hz) 

(b) Phase 

Figure 3.33 

ARX Model Frequency Response Comparison 
Blue=Actuator, Red=[7 6 1], Magenta=[l 1 14 1], Black=[l 1111] 
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Measure 


Result 


Average Output Magnitude 


0.1906 


Average Residual Magnitude 


0.0073 


Percentage of Average Output 


3.8% 


Average Output Difference Magnitude 


0.0025 


Mean Square Fit 


0.1736 



Table 3.6 

Numerical Results for Actuator #2 



7. Model Structure Acceptance for Actuator Number Three 

The plots resulting from the acceptance process are listed in appendix B 
(Figures B.l through B.7). Based on a delay of nk= 1 from table 3.1 and resulting 
parameter number plot (see appendix B Figure B.l) the structures selected for evaluation 
are listed in table 3.7. 



Number of Parameters 


Structure [na nb nk] 


Notes 


25 


15 10 1 


Automatic 


24 


.15 9 1 


ADC Based 


24 


14 10 1 


Manual 


17 


8 9 1 


MDL Based 


15 


87 1 


Automatic 


13 


8 5 1 


Automatic 


13 


76 1 


Manual 



Table 3.7 

Evaluated Model Structures for Actuator #3 

A comparison of the selected structure, [7 6 1], and the next best models is shown in 
Figure 3.34. Of particular interest, the best structures were those resulting from non- 
filtered data. A possible explanation is the absence of high frequency disturbances which 
filtering is designed to remove. Also of note is the fact that during excitation of actuator 
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number three the noise emitted was distinctly different from the other actuators when 
excited. This could be the result of slight manufacturing differences. 




io 3 



Frequency (Hz) 
(a) Magnitude 




io 3 



Frequency (Hz) 

(b) Phase 

Figure 3.34 

ARX Model Frequency Response Comparison 
Blue=Actuator#3, Red=[7 6 1], Magenta=[15 15 1], Black=[7 13 1] 
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Measure 


Result 


Average Output Magnitude 


0.2298 


Average Residual Magnitude 


0.0613 


Percentage of Average Output 


26% 


Average Output Difference Magnitude 


0.0037 


Mean Square Fit 


0.2061 



Table 3.8 

Numerical Results For Actuator #3 



8. Model Structure Acceptance for Actuator Number Four 

The plots resulting from the acceptance process are listed in appendix B 
(Figures C.l through C.7). Based on a delay of nk= 1 from table 3.1 and resulting 
parameter number plot (see appendix B Figure C.l) the structures selected for evaluation 
are listed in table 3.9. 



Number of Parameters 


Structure [na nb nk ] 


Notes 


26 


11 152 


ADC Based 


20 


8 12 2 


MDL Based 


18 


10 8 2 


Automatic 


12 


4 8 2 


Automatic 


26 


11 15 1 


ADC Based 


24 


11 13 1 


MDL Based 


17 


89 1 


Automatic 


13 


76 1 


Manual 



Table 3.9 

Evaluated Model Structures 

Similar to the case of actuator number one, the algorithm selected a delay of two. 
However, the best results were achieved by using a lower delay. The penalty paid was 
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degraded auto and cross correlation functions. This was again taken as a design trade off. 
The ARX structure \na=l nb=6 nk= 1] again proved to be the best choice (see Figure 
3.35). Numerical results are contained in table 3 10. 




Frequency (Hz) 



(a) Magnitude 




Frequency (Hz) 

(b) Phase 

Figure 3.35 

ARX Model Frequency Response Comparison 
Blue=Actuator#4, Red=[7 6 1], Magenta=[l 1 13 1], Black-[1 1111] 
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Measure 


Result 


Average Output Magnitude 


0.3052 


Average Residual Magnitude 


0.0066 


Percentage of Average Output 


2.1% 


Average Output Difference Magnitude 


0.0036 


Mean Square Fit 


0.2662 



Table 3.10 

Numerical Results for Actuator #4 

9. Model Structure Acceptance for Actuator Number Five 

The plots resulting from the acceptance process are listed in appendix B 
(Figures D.l through D.7). Based on a delay of nk=\ from table 3.1 and resulting 
parameter number plot (see appendix B Figure D.l) the structures selected for evaluation 
are listed in table 3.12. 



Number of Parameters 


Structure [na nb nk ] 


Notes 


28 


13 15 2 


ADC Based 


27 


12 152 


MDL Based 


22 


10 12 2 


Automatic 


18 


1262 


Automatic 


28 


13 15 1 


ADC/MDL Based 


24 


159 1 


Automatic 


15 


87 1 


Automatic 


13 


76 1 


Manual 



Table 3. 1 1 

Evaluated Model Structures 

The validation of actuator number five is virtually the same as actuator number four (see 
appendix B Figures C.2 through C.7, D.2 through D.7 and tables 3.10 and 3.12). 
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Strut #5 Magnitude 




10 3 

Frequency (Hz) 



(a) Magnitude 




10 3 

Frequency (Hz) 

(b) Phase 

Figure 3.36 

ARX Model Frequency Response Comparison 
Blue=Actuator #5, Red=[7 6 1], Magenta=[15 9 1], Black=[13 15 lj 



51 



Measure 


Result 


Average Output Magnitude 


0.3070 


Average Residual Magnitude 


0.0078 


Percentage of Average Output 


2.5% 


Average Output Difference Magnitude 


0.0032 


Mean Square Fit 


0.2830 



Table 3.12 

Numerical Results for Actuator #5 

10. Model Structure Acceptance for Actuator Number Six 

The plots resulting from the acceptance process are listed in appendix B 
(Figures E.l through E.7). Based on a delay of nk=\ from table 3.1 and resulting 
parameter number plot (see appendix B Figure E.l) the structures selected for evaluation 
are listed in table 3.13. 



Number of Parameters 


Structure [na nb nk ] 


Notes 


30 


15 15 1 


ADC/MDL Based 


24 


12 12 1 


Automatic 


20 


119 1 


Automatic 


18 


99 1 


Automatic 


13 


76 1 


Manual 



Table 3.13 

Evaluated Model Structures 

The comparison of the best model structures is shown in Figure 3.37 . The [7 6 1] 
structure selected for the other five actuators is the best choice for actuator number six. 
Numerical results are listed in table 3.14. 
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Strut #6 Magnitude 




Frequency (Hz) 
(a) Magnitude 



Strut #6 Phase 




Frequency (Hz) 

(b)Phase 

Figure 3.37 

ARX Model Frequency Response Comparison 
Blue= Actuator #6, Red=[7 6 1], Magenta=[15 15 1], Black=[12 12 1] 
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Measure 


Result 


Average Output Magnitude 


0.1551 


Average Residual Magnitude 


0.0055 


Percentage of Average Output 


3.5% 


Average Output Difference Magnitude 


0.0023 


Mean Square Fit 


0.1316 



Table 3.14 

Numerical Results for Actuator #6 
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D. APPLICATION OF SELECTED STRUCTURE 



The selection of the same basic model structure na-1 nb=6 nk= 1 for all six 
actuators is an indication that the plant dynamics of the actuators can be assumed to be 
identical. With the exception of actuator number three, all received the best results from 
prefiltered data. The design trade-offs were the degraded auto and cross correlation 
functions and the relative proximity of two poles and two zeros. This basic structure was 
used to develop more accurate ARMAX and Box-Jenkins models. Although the Box- 
Jenkins is the most accurate parametric estimation method, it was anticipated that the 
ARMAX would give the best results based on the nature of disturbances and on when 
they enter into the system (see Figure 3.3). 

The application and evaluation process is the same as that used to evaluate 
selected structures on the ARX model. A second order noise model was used for both the 
ARMAX and Box-Jenkins models. The standard choices are 2 nd and 4 th order models 
[Ref. 11]. Both were applied and the 2 nd order yielded a better result. After the 
application of the ARMAX and Box-Jenkins based models the overall performance of the 
ARX and ARMAX models was decidedly better. The only exception was on actuator 
number two. This will be addressed in section two. Based on the results obtained thus 
far indicating the actuator plants are essentially the same and the excellent performance of 
the ARX and ARMAX models, only the evaluation plots for the ARMAX model are 
included and can be found in appendix C. Summary and comparison results are listed for 
each actuator. 



1. Actuator Number One 

A graphical comparison of ARX, ARMAX and Box-Jenkins models is 
shown in Figure 3.38. Validation plots for the ARMAX [7 6 2 1] model are included as 
Figures A.l through A. 6 in appendix C. Numerical results are contained in table 3.15. 
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Strut #1 Magnitude 




(a) Magnitude 




10 3 



Frequency (Hz) 

(b) Phase 

Figure 3.38 

Model Frequency Response Comparison 

Blue= Actuator #1, Red=ARMAX[7 6 2 1], Magenta=ARX[7 6 1], Black=BJ[6 2 2 7 1] 
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Measure 


ARX[7 6 1] 
(Data Prefiltered) 


ARMAX[7 6 2 1] 
(Data Prefiltered) 


BJ[6 2 2 7 1] 
(Data Prefiltered) 


Ave. Output Magnitude 


0.1870 


0.1870 


0.1870 


Ave. Residual Magnitude 


0.0064 


0.0036 


0.0094 


Percent of Output Mag. 


3.4% 


2% 


5% 


Ave. Output Difference 


0.0024 


0.002 


0.0025 


Mean Square Fit 


0.1562 


0.1589 


0.1568 



Table 3.15 

Numerical Results for Actuator #1 



Based upon the above results the ARX[7 6 1] model was determined to best fit the 
performance of actuator number one. The choice between the ARX and ARMAX is a 
close one. This indicates the assumption that the disturbance or noise entered the system 
early and is accurately modeled by the dynamics (poles) of the plant. The addition of the 
C(q,0) polynomial is not necessary. The improved performance of the ARX and 
ARMAX models over the Box -Jenkins further supports the assumption regarding the 
nature of the disturbance. 



2. Actuator Number Two 

A graphical comparison of ARX, ARMAX and Box-Jenkins models is 
shown in Figure 3.39. Validation plots for the ARMAX [7 6 2 1] model are included as 
Figures B.l through B.6 in appendix C. Numerical results are contained in table 3.16. 
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Strut #2 Magnitude 




10 3 

Frequency (Hz) 



(a) Magnitude 




10 3 



Frequency (Hz) 

(b) Phase 

Figure 3.39 

Model Frequency Response Comparison 

Blue= Actuator #2, Red=ARMAX[7 6 2 1], Magenta=ARX[7 6 1], Black=BJ[6 2 2 7 1] 
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Measure 


ARX[7 6 1] 
(Data Pre filtered) 


ARMAX[7 6 2 1] 
(Data Prefiltered) 


BJ[6 2 2 7 1] 
(Data Prefiltered) 


Ave. Output Magnitude 


0.1906 


0.1906 


0.1906 


Ave. Residual Magnitude 


0.0073 


0.003 


0.0114 


Percent of Output Mag. 


3.8% 


1.6% 


6% 


Ave. Output Difference 


0.0025 


0.0027 


0.0023 


Mean Square Fit 


0.1736 


0.1705 


0.1616 



Table 3.16 

Numerical Results for Actuator #2 

The Box-Jenkins model showed excellent frequency response correlation, auto and cross 
correlation functions and some improved numerical results. However, due to the close 
performance of the ARMAX and Box-Jenkins models, the large uncertainty ellipses in 
the zero-pole plot (see appendix C) the ARMAX model was selected. 



3. Actuator Number Three 

A graphical comparison of ARX, ARMAX and Box-Jenkins models is 
shown in Figure 3.40. Validation plots for the ARMAX [7 6 2 1] model are included as 
Figures C.l through C.6 in appendix C. Numerical results are contained in table 3.17. 
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Strut #3 Magnitude 




10 3 



Frequency (Hz) 
(a) Magnitude 




10 3 10 4 



Frequency (Hz) 

(b) Phase 

Figure 3.40 

Model Frequency Response Comparison 

Blue=Actuator #3, Red=ARMAX[7 6 2 1], Magenta=ARX[7 6 1], Black=BJ[6 2 2 7 1] 
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Measure 


ARX [7 6 1] 
(Data not filtered) 


ARMAX[7 6 2 1] 
(Data Prefiltered) 


BJ[6 2 2 7 1] 
(Data Prefiltered) 


Ave. Output Magnitude 


0.2298 


0.2298 


0.2298 


Ave. Residual Mag. 


0.0613 


0.0040 


0.0215 


Percent of Output Mag. 


26% 


1.7% 


9.3% 


Ave. Output Difference 


0.0037 


0.0040 


0.0036 


Mean Square Fit 


0.2061 


0.2031 


0.2000 



Table 3.17 

Numerical Results for Actuator #3 

Although the non-filtered ARX model was chosen over the prefiltered ARX model, it 
appears that the performance of the prefiltered ARMAX is much better. This can be 
explained by possible manufacturing variation in actuator number three which requires 
the use of the C(q,0) polynomial to model the disturbance. 



4. Actuator Number Four 

A graphical comparison of ARX, ARMAX and Box-Jenkins models is 
shown in Figure 3.41. Validation plots for the ARMAX [7 6 2 1] model are included as 
Figures D.l through D.6 in appendix C. Numerical results are contained in table 3.18. 
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Strut #4 Magnitude 




Frequency (Hz) 
(a) Magnitude 




10 3 



Frequency (Hz) 

(b) Phase 

Figure 3.41 

Model Frequency Response Comparison 

Blue= Actuator #4, Red=ARMAX[7 6 2 1], Magenta=ARX[7 6 1], Black=BJ[6 2 2 7 1] 
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Measure 


ARX[7 6 1] 
(Data Prefiltered) 


ARMAX[7 6 2 1] 
(Data Not Prefiltered) 


BJ[6 2 2 7 1] 
(Data Prefiltered) 


Ave. Output Magnitude 


0.3052 


0.3052 


0.3052 


Ave. Residual Mag. 


0.0066 


0.0335 


0.0099 


Percent of Output Mag. 


2.1% 


11% 


3.2% 


Ave. Output Difference 


0.0036 


0.0034 


0.0034 


Mean Square Fit 


0.2662 


0.2661 


0.2572 



Table 3.18 

Numerical Results for Actuator #4 

This case appears similar to actuator number one where the disturbance is accurately 
characterized by the system pole model without the need for the C(q,0) polynomial. 



5. Actuator Number Five 

A graphical comparison of ARX, ARMAX and Box-Jenkins models is 
shown in Figure 3.42. Validation plots for the ARMAX [7 6 2 1] model are included as 
Figures E. 1 through E.6 in appendix C. Numerical results are contained in table 3.19. 
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Stmt #5 Magnitude 




10 2 10 3 



Frequency (Hz) 
(a) Magnitude 




to 3 



Frequency (Hz) 

(b) Phase 

Figure 3.42 

Model Frequency Response Comparison 
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Measure 


ARX[7 6 1] 
(Data filtered) 


ARMAX[7 6 2 1] 
(Data Prefiltered) 


BJ[6 2 2 7 1] 
(Data Prefiltered) 


Ave. Output Magnitude 


0.3070 


0.3070 


0.3070 


Ave. Residual Mag. 


0.0078 


0.0035 


0.0151 


Percent of Output Mag. 


2.5% 


1.1% 


4.9% 


Ave. Output Difference 


0.0032 


0.0033 


0.0034 


Mean Square Fit 


0.2830 


0.2688 


0.2698 



Table 3.19 

Numerical Results for Actuator #5 

6. Actuator Number Six 

A graphical comparison of ARX, ARMAX and Box-Jenkins models is 
shown in Figure 3.43. Validation plots for the ARMAX [7 6 2 1] model are included as 
Figures F.l through F.6 in appendix C. Numerical results are contained in table 3.20. 
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Strut #6 Magnitude 




10 3 

Frequency (Hz) 
(a) Magnitude 




10 3 



Frequency (Hz) 

(b) Phase 

Figure 3.43 

Model Frequency Response Comparison 

Blue= Actuator #6, Red=ARMAX[7 6 2 1], Magenta=ARX[7 6 1], Black=BJ[6 2 2 7 1] 
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Measure 


ARX[7 6 1] 
(Data not filtered) 


ARMAX[7 6 2 1] 
(Data Prefiltered) 


BJ[6 2 2 7 1] 
(Data Prefiltered) 


Ave. Output Magnitude 


0.1551 


0.1551 


0.1551 


Ave. Residual Mag. 


0.0055 


0.0035 


0.0086 


Percent of Output Mag. 


3.5% 


2.2% 


5.5% 


Ave. Output Difference 


0.0023 


0.0020 


0.0019 


Mean Square Fit 


0.1316 


0.1312 


0.1356 



Table 3.20 

Numerical Results for Actuator #6 



E. SUMMARY OF RESULTS 

Table 3.21 presents a summary of the model selection process. 



Actuator 


i 


2 


3 


4 


5 


6 


Model 

[Structure] 


ARX 
[7 6 1] 


ARM AX 
[7 6 2 1] 


ARM AX 
[7 6 2 1] 


ARX 
[7 6 1] 


ARM AX 
[7 6 2 1] 


ARM AX 
[7 6 2 1] 


Ave. Out. Mag. 


0.1870 


0.1906 


0.2298 


0.3052 


0.3070 


0.1551 


Ave. Res. Mag. 


0.0064 


0.003 


0.0040 


0.0066 


0.0035 


0.0035 


% of Out. 
Mag. 


3.4% 


1.6% 


1.7% 


2.1% 


1.1% 


2.2% 


Ave. Out. 
Dive. 


0.0024 


0.0027 


0.0040 


0.0036 


0.0033 


0.0020 


Mean Sq. Fit 


0.1562 


0.1705 


0.2031 


0.2662 


0.2688 


0.1312 



Table 3.21 

Summary of Numerical Results 



F. POLES AND ZEROS 

The ultimate goal in identifying the plant of the platform is to obtain the poles and 
zeros of the plant transfer function which can be used in design of an active vibration 
control law for the platform. Tables 3.22 and 3.23 present the poles and zeros extracted 
from the parameter estimation algorithm listed in appendix A. 
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1 


2 


3 


0. 1252+0.9 199i 


0.1080+0.9667i 


0.1198+0.9204i 


0. 1252-0.9 199i 


0. 1080-0. 9667i 


0.1198-0.9204i 


1.0137 


1.0222 


1.0204 


0.5371+0.6167i 


0.5887+0.6486i 


0.4925+0.6206i 


0.5371 -0.6 167i 


0.5887-0.6486i 


0.4925-0.6206i 



(a) Actuators #1 through #3 



4 


5 


6 


0.0674 + 0.6235i 


0.1971+0.8593i 


0.228 l+0.7945i 


0.0674 + 0.6235i 


0.1971-0.8593i 


0.228 l-0.7945i 


-1.2095 


1.0098 


1.0032 


1.0010 


-0.5566 


-0.1319+0.6080i 


-0.5568 


0.2913 


-0.131 9-0. 6080i 



(b) Actuators #4 through #6 
Table 3.22 

Actuator Transfer Function Poles 
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1 


2 


3 


0.1661 +0.9362i 


0. 1 122 + 0.9803i 


0. 1410 + 0.9402i 


0.1661 -0.9362i 


0.1 122 - 0.9803i 


0.1410 + 0.9402i 


0.8743 


0.9035 


0.7231 


0.5932 + 0.7591i 


0.5645 + 0.8101i 


0.5584 +0.777 li 


0.5932 -0.759H 


0.5645 -0.8101i 


0.5584 -0.777 li 


0.7101 +0.5677i 


0.7220 + 0.6035i 


0.5562 + 0.5203i 


0.7101 +0.5677i 


0.7220 - 0.6035i 


0.5562 + 0.5203i 



(a) Actuators #1 through #3 



4 


5 


6 


0.6905 + 0.6666i 


0.1657 + 0.9000i 


0.1 138 + 0.8029i 


0.6905 - 0.6666i 


0.1657 -0.9000i 


0.1 138 - 0.8029i 


0.6386 


0.8267 


0.6405 


-0.3437 + 0.4468i 


0.6710 + 0.6870i 


0.5522 + 0.7348i 


-0.3437 - 0.4468i 


0.6710 -0.6870i 


0.5522 -0.7348i 


0.2640 + 0.453 li 


0.4532 + 0.4747i 


0.3908 + 0.4466i 


0.2640- 0.453 li 


0.4532 - 0.4747i 


0.3908 -0.4466i 



(b) Actuators #4 through #6 
Table 3.23 

Actuator Transfer Function Zeros 



A review of the poles and zeros shows that zeros are very similar, especially among 
actuators one, two, three and six. The magnitude of the real parts of zeros for actuators 
four and five appear to vary more than the others. This is also true for the poles. This 
was unexpected and no explanation is offered. Any expected difference would have been 
attributed to actuator number three based on previous results and analysis. However, 
actuator number three is in relative concurrence with the other actuators. 
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G. TRANSFER FUNCTIONS 



The actuator transfer functions are as follows; 

ACTUATOR Number One 

13.53625 6 -31.64935 s + 42.53 1 75 4 - 39.48865 3 + 22.80125 2 -7.90855 
5 7 -3.81325 6 + 7.77935 s - 10.35325 4 + 9.71 025 3 - 6.39075 2 + 2.73365 - 0.6064 



ACTUATOR Number Two 

1 5.1 7095 6 - 36.64 645 5 + 5 1.4596 5 4 - 49.92965 3 + 30.8599 5 2 - 1 1.25845 
5 7 - 3.70105 6 + 7.56925 s -10.251 15 4 + 9.94825 3 - 6.87625 2 + 3.13325 - 0.7594 



ACTUATOR Number Three 

2 1.68755 6 - 48.69045 s + 64.51 985 4 - 59.84855 3 + 33.83705 2 -1 1.96735 
5 7 - 3.23425 6 + 6.06795 s - 7.50225 4 + 6.63805 3 - 4.1 175 5 2 + 1.64545 - 0.3361 



ACTUATOR Number Four 

13.48755 6 -30.79935 s +40.40325 4 - 36.90175 3 + 21.22365 2 - 7.54895 
5 7 - 3.92705 6 + 8.09535 s -10.80765 4 + 10.1 3385 3 - 6.63485 2 + 2.80155 - 0.6022 



ACTUATOR Number Five 

8.5 13 15 6 + 3.98235 s -9.42375 4 - 4.22425 3 + 0.930025 2 - 0.19505 
5 7 -2.03595 6 + 1.7435 s -0.48245 4 +0.01 755 3 -0.27735 2 +0.25775-0.1089 



ACTUATOR Number Six 

9.67685 6 - 1 1.56925 s + 1 1.059 15 4 - 9.1 858 5 3 + 2.52275 2 - 2.567 15 

5 7 -2.75425 6 +4.501 15 s -4.77445 4 + 3.65825 3 -1.96925 2 + 0.681 15 - 0.1253 

There is greater concurrence among actuators one, two, three and six. However, overall 
comparison of zeros, poles and transfer functions indicates that all six plants are very 
similar. 



70 



IV. ACTUATOR COUPLING AND VIBRATION CONTROL 

APPLICATION 

A. ACTUATOR COUPLING 

With the dynamic models constructed in Chapter Three it is now possible to model 
and investigate the interaction between actuators to see if indeed they can be controlled 
independently The basic ARX method has been proved to be quite accurate in modeling 
the plant transfer function of each actuator To determine the presence of actuator 
coupling each actuator was excited with 50mV pink noise and the output response of all 
six actuators measured. The input-output data gathered was introduced to the ARX[7 6 
1] model and the resulting model frequency response was plotted. Figures 4 1 through 4.5 
show the results from the excitation of actuator number one. 
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Figure 4 1 

Frequency Response Actuator #1 Vs. Actuator #2 
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Degrees 



60 



Magintude: lnput=50mV Pink Noise 




.40 1 1 — I I l Mi l 1 — l l M m 1 — l - 1 I M ill 

10 1 10 2 10 3 10 4 

Frequency (Hz) 

(a) Magnitude 




Frequency (Hz) 

(b) Phase 

Figure 4.2 

Frequency Response Actuator #1 Vs. Actuator #3 
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Magintude: lnput=50mV Pink Noise 




Frequency (Hz) 



(a)Magnitude 




Frequency (Hz) 

(b) Phase 

Figure 4.3 

Frequency Response Actuator #1 Vs. Actuator #4 
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Magintude: lnput=50mV Pink Noise 




Frequency (Hz) 



(a) Magnitude 




Frequency (Hz) 

(b) Phase 

Figure 4.4 

Frequency Response Actuator #1 Vs. Actuator #5 
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Degrees 




Frequency (Hz) 



(a) Magnitude 




Frequency (Hz) 

(b) Phase 

Figure 4.5 

Frequency Response Actuator #1 Vs. Actuator #6 
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In Figures 4. 1 through 4.4 there is discernible correlation between the frequency 
response of actuator number one and the frequency response of actuators two through five 
to the excitation of actuator number one. Of particular interest is the response at the 
resonance frequency of 1 .4 kHz. On average there is approximately a 10 dB difference at 
1.4 kHz for these actuators. In the frequencies below 1kHz the response is in the 
negative dB range. This would indicate that the passive stage is performing well in 
eliminating the low frequency disturbances. Although coupling is present, it is relatively 
low in comparison to the response of actuator one. However this is not the case with 
actuator number six. As in actuators two through five the response below 1 k is small. 
Above 1kHz, particularly at 1.4 kHz the response is close to a mirror image of the 
response of actuator number one. This suggests that there is significant high frequency 
coupling 



B. CONTROL APPLICATION AND COUPLING VERIFICATION 

Given the results of the previous section it is necessary to verify on the UQP. 
This was accomplished by applying a control law to the platform and checking for any 
evidence of coupling. The first evaluation conducted consisted of controlling an 
individual actuator and the second applied control to all six actuators simultaneously. 
There are numerous control candidates for use in A VC. For random broad band 
noise/vibrations, feedback control is preferable. Feed forward controllers are good 
candidates when a specific frequency of the disturbance is known. Local (decoupled) 
force feedback [Ref. 5], adaptive vibration control using two least mean square filters 
[Ref. 4], and absolute velocity feedback [Ref. 6] are examples that have been used in 
AVC. The latter of the three also used a geophone for sensing. The control law currently 
programmed for the system is a lead-lag compensator. This type of compensator is a 
variant of derivative and integral control. Lead compensation has the effect of lowering 
rise time and decreasing overshoot while lag compensation is used mainly to gain steady 
state accuracy of the system[Ref. 7]. A lead-lag controller was provided by CSA 
Engineering and was used to ensure the operation of each individual stmt. This control 



76 



law was modified to apply the same lead-lag compensator to all six actuators 
simultaneously based on the initial assumption of independently controllable actuators. 
The results previously obtained indicate that this cannot be done. 

1 . Compensator T ransfer F unction 

A Bode Plot of the magnitude and phase of the compensator (See Figure 
4.6) shows an infinite gain margin and a phase margin of approximately 60 degrees. 




Figure 4.6 

Compensator Bode Plot 
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A bode plot of the open loop tranfer function is provided in Figure 4.7. 




Figure 4.7 

Open Loop Transfer Function (Compensator* Actuator) 

Finally, a magnitude plot in terms of dB of closed to open loop is shown in Figure 4.8. 
From this plot it is observed that the best performance for this compensator will be 
between 12 Hz to 48 Hz with maximum reduction of 32dB at 25 Hz. A reduction of 
approximately 25 dB can be seen at 42 Hz. This controller was designed with an 
approximate bandwidth of 100 Hz. The spike at approximately 10 Hz is a faulty value 
resulting from system noise or an error in data collection. 
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2. Control Application 

The controller described in section one was first applied to actuator 
number two and then the .modified controller was applied to all six actuators 
simultaneously. Table 4.1 details the equipment preparation and settings. 
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Equipment 


Setting 


HP 33120A Function Generator 


Waveform Sinusoid 

Frequency 42 Hz 

Amplitude 500mv peak- 

to-peak 


KEPCO Bipolar Operational Power 
Supply/Amplifier 


Current Control Mode +/- 1 Amp 



Table 4.1 

Equipment Settings 



3. Experiment 

The first task is to look at the platform in open loop mode. From this the 
local environment can be assessed. The plots in Figures 4.9 and 4.10 show what is 
sensed by the geophone on actuator two with the 42 Hz disturbance and no control (open 
loop). 



Open Loop-42Hz Disturbance 




Figure 4.9 

Ambient PSD (0-5kHz) 
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Open Loop-42Hz Disturbance 




Figure 4.10 

Ambient PSD (0-200Hz) 

Of note in Figure 4.9 is the 1.4 kHz resonance frequency. This is the result of power 
being applied to the actuator from the HVPS. At approximately 30 Hz there is an 
environmental disturbance of unknown origin. In the open loop mode the 42Hz 
disturbance is clearly visible (see Figure 4.10). Figure 4.11 shows the power spectral 
density (PSD) of the sensor output for actuator number two for open and closed loop 
operations (both single and all six actuators). With only power applied to the 
piezoceramic actuator (no control), excitation of the resonance frequency is visible as the 
blue trace in Figure 4.11(a). Closing the control loop for single and multi-actuator 
control increases the excitation. Figure 4.12 shows the sensor output for actuators two, 
three, and six with control applied to actuator number two. Actuator number three shares 
the same base node with actuator number two. The coupling between actuators two and 
three can be seen in Figure 4.12(a). The application of control excites the resonance of 
all three actuators, but the interaction between actuators two and three is the most 
pronounced. 
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Blue-Open Loop; Magenta-Strut #2 Closed Loop; Red-All 6 Struts Closed Loop 




(a) PSD (0-5kHz) Sensed By Geophone on Actuator #2 
Blue-Open Loop; Magenta-Strut #2 Closed Loop; Red-All 6 Struts Closed Loop 




Frequency 

(b) PSD (0-200Hz) Sensed By Geophone on Actuator #2 

Figure 4. 1 1 

Open Vs. Closed Loop Performance PSD 
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Blue-#2 Geophone; Magenta-#3 Geophone; Red-# 6 Geophone 




(a) PSD (0-5kHz) Sensor Outputs 
Blue-#2 Geophone; Magenta-#3 Geophone; Red-# 6 Geophone 




Frequency 



(b) PSD (0-200Hz) Sensor Outputs 
Figure 4.12 

Actuator #2 Closed Loop Performance PSD 
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V. SUMMARY 



A. CONCLUSIONS 

The process of system identification provided valuable insight into the UQP. By 
treating the UQP as a “Black Box” and using ready-made models the complex system 
dynamics were overcome and very accurate models of the system were developed. Some 
possible reasons for differences between the models and the actual system can be 
explained by non-removal of faulty values. Due to the number of samples used to build 
and validate the models, hand smoothing of the data was not practical. Based on the 
average residuals for each actuator being less than 4% of the average model output, 
outliers did not pose a significant problem. In comparing the validation data (numerical 
and graphical), actual frequency responses, and zeros, poles and transfer functions there 
are distinct similarities. The fact that all six are described by seven poles and six zeros 
was anticipated. However, some variation among actuators was expected. There is 
strong correlation of actuators one, two, three and six are. Although the numerical values 
for the zeros, poles and transfer function of actuators four and five did not correlate with 
the other actuators as strongly, they still retained the same pole-zero structure. Based on 
this it can be stated that for the purposes of characterizing the system and control design, 
all six actuators are identical. 

The sensitivity of the geophones resulting in the introduction of external 
disturbances could be a source of inaccuracy, however, they are considered small in 
comparison to the disturbances generated by the other actuators. The emergence of ARX 
and ARMAX models as the most accurate models validates the assumption that these 
disturbances to the system entered the process early on. The disturbance created by the 
excitation of an actuator in turn causes excitation of other actuators, especially for 
actuators that share the same aluminum base in the region above 200 Hz. Since the 
disturbance is caused by another actuator, the choice of modeling the noise with the 
system dynamics (poles) was an accurate one. The interaction or coupling among 
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actuators, particularly between base adjoining actuators is distinct. Although the passive 
damping mechanism works well in eliminating low frequency interaction, high frequency 
disturbances are readily transmitted between actuators. For actuators not sharing a 
common base node this interaction is minimal and could possibly be ignored. However, 
for base node sharing elements this problem is significant. 

The lead-lag compensator worked well for the single actuator in the 1Hz to 100Hz 
range, providing over 20dB damping of the 42Hz disturbance. However, above 200 Hz 
the controller exacerbates the disturbances. This problem is amplified in the case where 
all six actuators are controlled simultaneously. The interaction between actuators 
compounds the high frequency degraded performance of the controller. This causes it to 
go unstable. The high frequency interaction between base sharing nodes must be 
addressed before attempting to control all six actuators simultaneously. 

B. RECOMMENDATIONS FOR FUTURE WORK 

With accurate plant models of the actuators the next logical step is control law 
design. With knowledge of the system it is possible to design controllers using either 
classical or modem techniques. The main concern in the design is the high frequency 
interaction in the range of 1kHz to 1.5 kHz. The revelation of actuator interaction does 
not necessarily mean that it is impossible to control each actuator independently. The 
dynamic compensator applied to an individual actuator worked well below 100Hz but did 
not roll off sufficiently. This resulted in high frequency instability in control of all six 
actuators. With knowledge of the plant transfer function it is possible to design a 
dynamic compensator capable of controlling each actuator independently. However, the 
design must be such that it does not excite the 1 .4 kHz resonance. On the other hand, it 
may be necessary to treat the UQP as a multiple input/multiple output (MIMO) system 
and a modem control design is therefore required. Another possible remedy to the 
interaction between base node sharing elements is to change the material of the 
attachment node. The use of aluminum acts as an excellent conductor of vibration from 
one actuator to another. By using vibration absorbing material such as high stiffness 
rubber it may be possible to reduce the interaction sufficiently such that each actuator 
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could be controlled independently using single input/single output (SISO) classical design 
techniques. Replacing the upper node material could also reduce the overall transmission 
of unwanted vibration, not only from interaction but external sources as well. 
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APPENDIX A. COMPUTER PROGRAMS 



All computer codes included in this thesis were written by Naval Postgraduate 
School Space Research and Design Center unless otherwise noted. The code included in 
this appendix is organized along lines of the thesis. 

A. SYSTEM IDENTIFICATION AND MODEL CONSTRUCTION 
% sysident.m : MAIN PROGRAM 

% Written by: George D. Beavers June 1997 

% This MATLAB “.m” file reads in the user requested actuator input output data 
% from the “4h*.mat” files. “*” represents the actuator that is being excited by the 
% 50mV “pink noise”. It calls the user written functions 

% “arxmod.m”, “armaxmod.m” and “bjmod.m” to calculate models based on the 
% ARX, ARMAX and Box-Jenkins parameter estimation methods respectively 
% The zeros, poles and transfer functions are then extracted from the models 
% returned by the user written functions. 

% Each “4h*.mat” file has 7 rows of data 

% rows 1-6 = geophone output from the respective actuator (rowl=actuator #1) 

% row 7 is the 50mv “pink noise” input 

% Information on the MATLAB SITB functions can be found in reference[12] 

clear 

dum='y’; 

while (dum=='y'), 
j=input('Enter strut number » '); 
disp(['load 4h' int2str(j)]) 

eval(['load ' '4h' int2str(j) ]) %Load input output data % 

z2=[trace_y(j, 1:10000)' trace_y(7,l: 10000)']; %Pull out requested channel % 

z2=dtrend(z2); %Remove trends & adjust levels% 

z3=[trace_y(j,l 0001 :20000)' trace_y(7, 1000 1 :20000)'] ; 
z3=dtrend(z3); 

disp(['Choose Structure Selection Method']) 
disp([' 1= Automatic']) 
disp(['2=Manual’]) 
temp=input('»> '); 

disp(['Choose a Parametric Estimation Model']) 
disp(['arx =1'; 'armax =2'; 'bj =3']) 
chc=input('»> ') 

% The following subroutine is used to determine input delay and number of parameters'^ 
% and structure na,nb based on that delay for an ARX based model% 
if temp==l 

V=arxstruct(z2,z3,struc(2,2, 1 :5)); 



% Use algorithm to determine % 
%[na, nb, nk] or input them % 
% manually % 
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[nn, V m]=selstruc(V,' AIC) 
inp=input(’enter delay') 

V=arxstruct(z2,z3,struc(l : 15, 1 : 15,inp)); 

nnA=selstruct(V,'AIC') 

nnM=selstruct(V,'MDL') 

Figure(l) 

nns=selstruct(V) % Generate parameter Vs. Loss fcn% 

% plot % 

% Determine which model to construct then call the function% 

% Using the automatically generated ARX based Structure % 
if chc==l 

arxmod(j,nns,z2,z3); 
elseif chc==2 

nc=input('Enter Order of Noise model nc »'); 
itr=input('Enter Maximum Number of Iterations »'); 
narm=[nns(l,l) nns(l,2) nc nns(l,3)]; 
armaxmod(j,narm,itr,z2,z3) 

else 

on=input('Enter Order of Noise model [nc nd] »') 
itr=input('Enter Maximum Number of Iterations »') 
nbj=[nns(l,2) on nns(l,l) nns(l,3)] 
bjmod(j,nbj,itr,z2,z3) 

end 

% Generate models based on manually entered structure % 
else 



if chc==l 

nns=input('Enter Structure [na nb nk] »'); 

[arxm,arxmft]=arxmod 1 (j,nns,z2,z3); 

[zepo 1 ,k 1 ]=th2zp(arx m) % Extracts zeros-poles 

[num 1 ,den 1 ]=th2tf(arxm) % Extracts transfer function 

[zepolf,klf]=th2zp(arxmft) % Extracts zeros-poles (filtered data) 

[num 1 f,den 1 f]=th2tf(arxmft) % Extracts tranfer function(filtered data) 
present(arxm) % Extracts parameters 

present(arxmft) % Extracts parameters(filtered data) 

elseif chc=2 

narm=input('Enter Structure [na nb nc nk] »'); 
itr=input('Enter Maximum Number of Iterations »') 
[armxm,armxmft]=armaxmdf(j,narm,itr,z2,z3); 
[zepo2,k2]=th2zp(armxm) % Extracts zeros-poles 

[num2,den2]=th2tf(armxm) % Extracts tranfer function 

[zepo2f,k2f]=th2zp(armxmft) % Extracts zeros-poles (filtered data) 
[num2f,den2f]=th2tf(armxmft) % Extracts tranfer function(filtered data) 
present(armxm) % Extracts parameters 
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present(armxmft) % Extracts parameters(filtered data) 

else 

nbj=input('Enter Structure [nb nc nd nf nk] »') 
itr=input('Enter Maximum Number of Iterations »') 
[bjm,bjmft]=bjmdf(j,nbj,itr,z2,z3); 

end 

end 

dum=input('Continue ?» Vs'); 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



% arxmod.m % 

% ARX Model Construction and Analysis % 

% George D. Beavers % 

% June 1997 % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% This program generates a model based on the input-output data of a system using % 
% the ARX parameter estimation method. It also generated a series of graphical and % 
% numerical results in order to evaluate the suitability of the model % 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

function [arxm,arxmft]=arxmod(j,nns,z2,z3) 

%%%%%%%%%%%%%%%%%%%%%%%%%% 



% arxm=ARX based model % 

% arxmft=ARX based model (filtered data) % 

% j=actuator % 

% nns=structure([na nb nc]) % 

% z2=input-output data for model construction % 

% Z3=input-ouput data for model evaluation % 



%%%%%%%%%%%%%%%%%%%%%%%%%% 

na=nns(l,l); 

NA=int2str(na); 

nb=nns(l,2); 

NB=int2str(nb); 

nk=nns(l,3); 

NK=int2str(nk); 

<^****4:************^,;^^ MODEL *************************** 

arxm=arx(z2,nns); 

arxm=sett(arxm,0.1e-3); % Set sample interval % 

garxm=th2ff(arxm); % Convert to frequency function% 
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% This section loads separate data to compute the actual actuator frequency response% 
% and then plot it against the model’s frequency response % 

Figure(2) 



eval(['load s',int2str(j),'m']) % 

Xl=o2ilx; 

Yl=20.*logl0(abs(o2il)); 

semilogx(X 1 , Y 1 ,'b'), 
hold on 

semilogx(garxm(:,l)/(2*pi),20*logl0(garxm(:,2)),'r'); 

title(['Blue=Actuator#',num2str(j),'Red=ARX MODEL[',NA,' ',NB,' ',NK,']']) 

xlabel('Frequency (Hz)') 

ylabel(’Magnitude(dB)') 

grid 

hold off 
Figure(3) 

eval(['load s',int2str(j),'p']) 

X2=o2ilx; 

Y2= 1 80*unwrap(angle(o2i l))/pi; 

semilogx(X2 , Y2,’b') 
hold on 

semilogx(garxm( : , 1 )/(2*pi),(garxm(: ,4)),'r') ; 

title(['Blue=Actuator#’mum2str(j),', Red=ARX MODEL[’,NA,' ',NB,' ',NK,']’]) 
xlabel('Frequency (Hz)') 
ylabel('Phase (Deg)') 
grid 

hold off 

<y 0 ******************fijfj£ MODEL (DATA PREF1LTERED)*****************% 
% The above process is repeated using prefdtered data. A 10 th order Butterworth filter % 
% is used. The passband is given in terms of fractions of the Nyquist frequency (0.52)% 
% This is a low pass filter with an upper bound of approximately 1.5kHz 

% 

zf=idfilt(z2,10,.52); %Prefilter model construction data% 

arxmft=arx(zf,nns); 
arxmft=sett(arxmft,0. le-3); 
garxmft=th2ff(arxmft); 

eval(['load s',int2str(j),'m']) 

Xl=o2ilx; 
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Y 1 =20. *log 1 0(abs(o2i 1 )); 



Figure(4) 

semilogx(X 1 ,Y 1 ,'b') 
hold on 

semilogx(garxmft( : , 1 )/(2*pi),20*log 1 0(garxmft(:,2)),’r’); 

title(['Blue=Actuator#',num2str(j),', Red=ARX MODEL[\NA,' ',NB,' ',NK,’] *Data 

Prefiltered*’]) 

xlabel('Frequency (Hz)') 

ylabel(’Magnitude(dB )') 

grid 

hold off 
Figure(5) 

eval(['load s’,int2str(j),'p']) 

X2=o2ilx; 

Y2= 1 80*unwrap(angle(o2i 1 ))/pi; 
semilogx(X2,Y2,'b') 

hold on 

serrulogx(garxrrTt(:,l)/(2*pi),(garxmft(:,4)),'r’); 

title(['Blue=Actuator#',num2str(j),', Red=ARX MODEL[’,NA,' ',NB,' ',NK,'] *Data 

Prefiltered*']) 

xlabel('Frequency (Hz)’) 

ylabel('Phase (Deg)') 

grid 

hold off 

% ****************** GRAPHICAL AND numerical analysis***********% 

% This section generates plots and numerical results used in evaluating the model % 
Figure(6) 

% Apply the input row of the validation data to the model and then compare with % 
% the actual output of the system corresponding to the validation data and plot results% 
[y h ,fit] =compare(z3,arxm) ; 
axis([5000,5 100,-1.0,1.0]) 

title(['Plot of ARX MODEL [',NA,' ',NB,' ’,NK,'] Output vs. Actual output']) 
text(5020,0.75,'Red=Model Output, Blue= Measured Output') 
xlabel(’S ample') 
ylabel('Output Magnitude') 

Figure(7) 

[yhf,fitfj=compare(zfv,arxmft); %Prefiltered data% 

title(['Plot of ARX MODEL [',NA,' ’,NB,’ ',NK,'] Output vs. Actual output *Data 

Prefiltered*']) 
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axis([5000,5 100,- 1 .0, 1 .0]) 

text(5020,0.75,'Red=Model Output, Blue= Measured Output') 

xlabel('Sample') 

ylabel('Output Magnitude’) 

% Calculate the residuals and plot the auto and cross correlation functions % 
Figure(8) 

el=resid(z3,arxm); % Non-filtered data 

Figure(9) 

e2=resid(zfv,arxmft); % Prefiltered data 

% Plot the residuals % 

Figure(lO) 

plot(el,'m'),grid 

title(['Residual of ARX Model [',NA,' ',NB,' ',NK,']’]) 

ylabel('Residual') 

xlabel('Sample') 

axis([5000,5 100,-0. 1,0.1]) 

Figure(ll) 

plot(e2,'m'),grid 

titl e( ['Residual of ARX Model [',NA,' ',NB,' ',NK,'] *Data Prefiltered*']) 

axis([5000,5 100,-0.05,0.05]) 

xlabel('Sample') 

ylabel('Residual') 



% Calculate Numerical Results% 
elpos=abs(el); 
e2pos=abs(e2); 
ypos=abs(z3(:,l)); 
ave=mean(e 1 pos) 
aveft=mean(e2pos) 
ou ta ve=mean (y pos) 
fit 
fitf 



% Average residual 
% Average residual 
% Average actual output 
%Mean Square Fit calculated above 
%Mean Square Fit calculated above 



% Calculate and plot the actual output-model output 

z3out=z3(:,l); 

for i=5000:5100; 

difft(i)=z3out(i)-yhf(i); 

end 

z3out=z3(:,l); 
for j=5000:5100; 

diff(j)=z3out(j)-yh(j); 

end 

Figure(14) 

plot(diff,'r') 
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axis([5000,5 100,-0.5,0.5]) 

titled Actual Output - ARX Model [',NA,' ’,NB,' ',NK,'] Model’]) 
xlabel(’Sample') 

Figure(15) 

plot(difft,'r') 

axis([5000,5 100,-0.5,0.5]) 

title([' Actual Output - ARX Model [',NA,' ,NB,' ',NK,'] Model *Data Prefiltered*']) 
xlabel('Sample') 

avediff=mean(abs(diff)) %Ave. difference between actual and model outputs 
avedft=mean(abs(difft)) % Ave. difference between actual and model outputs (filter) 

% Generate Zero-Pole plots 
Figure(12) 

zpplotl(th2zp(arxm),3,[],1.2) 

title(['Zero Pole Plot ARX Model [',NA/ ',NB,' ',NK,'] Model']) 

Figure(13) 

zpplot 1 (th2zp(arxmft) ,3, [] , 1 .2) 

title([’Zero Pole Plot ARX Model [',NA,' ',NB,' ',NK,'] Model *Data Prefiltered*']) 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



% armaxmod.m % 

% ARMAX Model Construction and Analysis % 

% George D. Beavers % 

% June 1997 % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% This program generates a model based on the input-output data of a system using % 
% the ARMAX parameter estimation method. It also generated a series of graphical % 
% and numerical results in order to evaluate the suitability of the model % 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

function [arxm,arxmft]=armaxmod(j,nns,itr,z2,z3) 

%%%%%%%%%%%%%%%%%%%%%%%%%% 



% arxm=ARMAX based model % 

% arxmft=ARMAX based model (filtered data) % 

% j=actuator % 

% nns=structure([na nb nc]) % 

% z2=input-output data for model construction % 

% Z3=input-ouput data for model evaluation % 



%%%%%%%%%%%%%%%%%%%%%%%%%% 

na=nns(l,l); 

NA=int2str(na); 

nb=nns(l,2); 

NB=int2str(nb); 
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nc=nns(l,3); 

NC=int2str(nc); 

nk=nns(l,4); 

NK=int2str(nk); 



<^*****************^pj^j^j^ MODEL **********************************% 

% This section loads separate data to compute the actual actuator frequency response% 
% and then plot it against the model’s frequency response % 

arxm=armax(z2,nns,itr); 
arxm=sett(arxm,0. le-3); 
garxm=th2ff (arxm) ; 



eval(['load s',int2str(j),'m']) 

Xl=o2ilx; 

Y 1 =20. *log 1 0(abs(o2i 1 )) ; 

Figure(2) 

semilogx(X 1 , Y 1 ,’b’), 
hold on 

semilogx(garxm( : , 1 )/(2*pi),20*log 1 0(garxm( : ,2)),'r') ; 

title(['Blue=Actuator#',num2str(j),'Red=ARMAX MODEL[',NA,' ’,NB,' ',NC,' ',NK,’]']) 

xlabel(’Frequency (Hz)') 

ylabel('Magnitude(dB)') 

grid 

hold off 

eval(['load s',int2str(j),'p’]) % Load actual actuator data (phase)% 

X2=o2ilx; 

Y2= 1 80*unwrap(angle(o2i 1 ))/pi; 

Figure(3) 

semilogx(X2,Y2,’b') 
hold on 

semilogx(garxm(: , 1 )/ (2*pi),(garxm(:,4)),'r'); 

title(['Blue=Actuator#',num2str(j),'Red=ARMAX MODEL[',NA,' ',NB,' \NC,' ’,NK,T]) 
xlabel('Frequency (Hz)’) 
ylabel(’Phase (Deg)') 
grid 

hold off 

<^*************^jyy| i ^^ MODEL (DATA PREFILTER£D********************% 
% The above process is repeated using prefiltered data. A 10 th order Butterworth filter % 
% is used. The passband is given in terms of fractions of the Nyquist frequency (0.52)% 



% Set sampling interval % 

% Convert to frequency function% 



% Load actual actuator data (magnitude)% 
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°Io This is a low pass filter with an upper bound of approximately 1 ,5kHz 
% 

zf=idfilt(z2,10,.52); % Prefilter model construction data % 

arxmft=armax(zf,nns,itr); 

arxmft=sett(arxmft,0, le-3); 

garxmft=th2ff(arxmft); 

eval(['load s',int2str(j),’m']) 

Xl=o2ilx; 

Y 1 =20. Hog 1 0(abs(o2i 1 )) ; 

Figure(4) 

semilogx(X 1 , Y 1 ,'m') 
hold on 

serrulogx(garxmft(:,l)/(2*pi),20*logl0(garxmft(:,2)),'g'); 

title(['Blue=Actuator#’,num2str(j),’Red=ARMAX MODEL[',NA,' ',NB,' ',NC,' ',NK,'] 

*Data Prefiltered*']) 

xlabel('Frequency (Hz)') 

ylabel('Magnitude(dB)') 

grid 

hold off 

eval(['load s',int2str(j),’p’]) 

X2=o2ilx; 

Y2=180*unwrap(angle(o2i l))/pi; 

Figure(5) 

semilogx(X2,Y2,'b') 
hold on 

semilogx(garxmft(:,l)/(2*pi),(garxmft(:,4)),’r’); 

title(['Blue=Actuator#',num2str(j),'Red=ARMAX MODEL[',NA,’ ',NB,' \NC,' ’,NK,'] 

*Data Prefiltered*']) 

xlabel('Frequency (Hz)') 

ylabel(’Phase (Deg)') 

grid 

hold off 

NUMERICAL ANALYSIS********** *% 
% This section generates plots and numerical results used in evaluating the model % 

% Apply the input row of the validation data to the model and then compare with % 

% the actual output of the system corresponding to the validation data and plot results 
% 

Figure(6) 

[yh,fit]=compare ] (z3,arxm); 
axis([5000,5 1 00,- 1 .0, 1 .0]) 
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title(['Plot of ARMAX MODEL[',NA; ',NB,' ’,NC,' ',NK,'] Output vs. Actual output’]) 

text(5020,0.75,’Red:Model Output, Green: Measured Output') 

grid 

Figure(7) 

[yhf,fitf]=compare 1 (zfv,arxmft); 

title(['Plot of ARMAX MODEL[',NA,' ',NB,' ’,NC,’ ',NK,'] Output vs. Actual output 

*Data Prefiltered*']) 

axis([5000,5 1 00,- 1 .0, 1 .0]) 

xlabel('Sample’) 

grid 

% Calculate residuals and generate and plot auto & cross correlation functions % 
Figure(8) 

el=residl(z3,arxm); 

Figure(9) 

e2=resid 1 (zfv,arxmft); 



% Plot Residuals % 

Figure(lO) 

plot(el,'m'),grid 

title( ['Residual of ARMAX MODEL[',NA,' ',NB,' ',NC,' ',NK,']']) 

ylabel('Residual') 

axis([5000, 5 100,-0. 1,0.1]) 

Figure(l 1) 
plot(e2,'m'),grid 

title(['Residual of ARMAX MODEL[',NA,' ',NB,' ’,NC,’ ',NK,'] *Data Prefiltered*’]) 
axis([5000,5 100,-0.05,0.05]) 
xlabel ('Sample') 
ylabel(’Residual') 



% Calculate Numerical Results% 
elpos=abs(el); 
e2pos=abs(e2); 
ypos=abs(z3(:,l)); 

ave=mean(elpos) % Average residual 

aveft=mean(e2pos) % Average residual (Prefiltered) 

outave=mean(ypos) % Average actual output 



fit 

fitf 



% Mean Square Fit calculated above 
% Mean Square Fit (Prefiltered) 



% Calculate actual output minus model output an plot results % 

z3out=z3(:,l); 
for i=5000:5100; 
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difft(i)=z3out(i)-yhf(i); 

end 

z3out=z3(:,l); 
for j=5000:5100; 

diff(j) = z3out(j)-yh(j); 

end 

Figure(14) 

plot(diff,’r) 

axis([5000, 5 100,-0.5,0.5]) 

title(['Actual Output - ARMAX[',NA,' ’,NB,' ’,NC,' ',NK,'] Model Output']) 
grid 

Figure(15) 

plot(difft,’r') 

axis([5000,5 100,-0.5,0.5]) 

title(['Actual Output - ARMAXf'^NA,' ',NB,' ’,NC,' ',NK,’]Model Output *Data 

Filtered*']) 

grid 

avediff=mean(abs(diff)) % Ave. output difference 

avedft=mean(abs(difft)) % Ave. output difference (Prefiltered) 

% Generate Zero-Pole Plots 
Figure(12) 

zpplot(th2zp(arxm),3,[],1.2) 

title(['Zero Pole Plot ARMAX MODEL[',NA,' ’,NB,' ',NC,' ’,NK,T]) 

Figure(13) 

zpplot(th2zp(arxmft) ,3 , [] ,1.2) 

title(['Zero Pole Plot ARMAX [',NA,' ',NB,' ’,NC,' ',NK,'] Model *Data Prefiltered*']) 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% bjmod.m % 

% Box-Jenkins Model Construction and Analysis % 

% George D. Beavers % 

% June 1997 % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% This program generates a model based on the input-output data of a system using % 
% the Box-Jenkins parameter estimation method. It also generated a series of graphical 
% and numerical results in order to evaluate the suitability of the model % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [bjm,bjmft]=bjmod(j,nns,itr,z2,z3) 
%%%%%%%%%%%%%%%%%%%%%%%%%% 



% bjm=B ox- Jenkins based model % 

% bjmft=Box -Jenkins based model (filtered data) % 

% j=actuator % 

% nns=structure([na nb nc]) % 

% itr=number of iterations % 
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% z2=input-output data for model construction % 

% Z3=input-ouput data for model evaluation % 

%%%%%%%%%%%%%%%%%%%%%%%%%% 
nb=nns(l,l); 

NB=int2str(nb); 

nc=nns(l,2); 

NC=int2str(nc); 

nd=nns(l,3); 

ND=int2str(nd); 

nf=nns(l,4); 

NF=int2str(nf); 

nk=nns(l,5); 

NK=int2str(nk); 



<^ *********** ******gQ^_j£jSjjQj s jg MODEL************************** *****% 
% This section loads separate data to compute the actual actuator frequency response % 
% and then plot it against the model’s frequency response % 



arxm=bj(z2,nns,itr); 
arxm=sett(arxm,0.1e-3); 
garxm=th2ff(arxm); 

eval(['load s',int2str(j),'m']) 

Xl=o2ilx; 

Yl=20.*logl0(abs(o2il)); 

Figure(2) 

semilogx(X 1 , Y 1 ,'b'), 
hold on 

senrulogxCgarxmO.iy^pi^O^loglOCgarxmO^Xh'r'); 

title(['Blue=Actuator#',num2str(j),'Red=Box-Jenkins MODEL[',NB,' ',NC,' ',ND,' ',NF,' 
’,NK,T]) 

xlabel('Frequency (Hz)’) 
ylabel('Magnitude(dB)') 
grid 

hold off 



%Build the model % 

% Set the sampling interval % 

% Convert to frequency function % 

%Load data for actuator frequency function % 



eval(['load s',int2str(j),'p'J) 

X2=o2ilx; 

Y2=180*unwrap(angle(o2il))/pi; 

Figure(3) 

semilogx(X2, Y 2,'b') 
hold on 

semilogx(garxm(:, 1 )/(2*pi),(garxm(:,4)),'r’); 

tit!e([’Blue=Actuator#',num2str(j),'Red=Box-Jenkins MODEL[’,NB,’ ',NC,' ',ND,' ',NF,' 
',NK,']']) 
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xlabel(’Frequency (Hz)') 
y label ('Phase (Deg)') 
grid 

hold off 

MODEL (DATA PREFILTERED******** ****** (^ 
% The above process is repeated using prefiltered data. A 10 th order Butterworth filter % 
% is used. The passband is given in terms of fractions of the Nyquist frequency (0.52)% 
% This is a low pass filter with an upper bound of approximately 1 .5kHz % 

zf=idfilt(z2, 10..52); 
arxmft=bj(zf,nns,itr); 
arxmft=sett(arxmft,0. 1 e-3); 
garxmft=th2ff(arxmft); 

eval(['load s',int2str(j),'m']) 

Xl=o2ilx; 

Y 1 =20. *log 1 0(abs(o2i 1 )) ; 

Figure(4) 

semilogx(X 1 , Y 1 ,'b') 
hold on 

semilogx(garxmft(:,l)/(2*pi),20*logl0(garxmft(:,2)),’r'); 

title(['Blue=Actuator#’,num2str(j),'Red=Box -Jenkins MODEL[',NB,' ',NC,' ',ND,' ',NF,' 
',NK,']*Data Prefiltered*']) 
xlabel('Frequency (Hz)') 
ylabel('Magnitude(dB)') 
grid 

hold off 

eval(['load s’,int2str(j),'p']) 

X2=o2ilx; 

Y2=180*unwrap(angle(o2il))/pi; 

Figure(5) 

semilogx(X2,Y2,'m') 
hold on 

semilogx(garxmft(:,l)/(2*pi),(garxmft(:,4)),'g'); 

title(['Blue=Actuator#',num2str(j),'Red=Box -Jenkins MODEL[',NB,' ’,NC,' ’,ND,’ ',NF,' 
',NK,']*Data Prefiltered*']) 
xlabel('Frequency (Hz)') 
ylabel('Phase (Deg)') 
grid 

hold off 



% Prefilter model construction data % 

% Construct model % 

% Set sampling interval % 

% Convert to frequency function % 

% Load data for actuator frequency response 
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% ****************** GI ^ HICAL AND NUMERICAL ANALYSIS***********% 
% This section generates plots and numerical results used in evaluating the model % 
% Comparison of model output to actual output % 

% Apply the input row of the validation data to the model and then compare with % 
% the actual output of the system corresponding to the validation data and plot results 
% 



Figure(6) 

[yh,fit]=compare(z3 ,arxm) ; 
axis([5000,5 100,-1.0,1.0]) 

title(['Plot of Box- Jenkins [',NB,' \NC,' \ND,’ ',NF,' \NK,'] Output vs. Actual output']) 

text(5020,0.75,'Red:Model Output, Blue: Measured Output') 

xlabel(’Sample') 

Figure(7) 

[yhf,fitf]=compare(z3,arxmft); 

title([’Plot of Box- Jenkins MODEL [',NB,' ’,NC,’ ',ND,' ',NF,' ’,NK,'] Output vs. Actual 
output *Data Prefiltered*']) 

text(5020,0.75,'Red:Model Output, Blue: Measured Output’) 

axis([5000,5 100,- 1 .0, 1 .0]) 

xlabel('Sample') 

% Calculate residuals and auto and cross correlation functions and plot results 
Figure(8) 

el=residl(zfv,arxm); 

Figure(9) 

e2=resid 1 (zfv,arxmft); 

Figure(lO) 

plot(el,'m'),grid 

title( ['Residual of Box-Jenkins Model [',NB,' ’,NC,’ ',ND,' ',NF,' ’,NK,’]']) 

ylabel('Residual') 

xlabel('Sample') 

axis([5000,5 100,-0. 1 ,0. 1 ]) 

Figure(l 1) 
plot(e2,'m'),grid 

title(['Residual of Box-Jenkins Model [',NB,' ’,NC,’ ',ND,' ’,NF,’ ',NK,'] *Data 
Prefiltered*']) 

axis([5000,5 100,-0.05,0.05]) 

xlabel(’Sample') 

ylabel(’Residual') 
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% Numerical results 

elpos=abs(el); 

e2pos=abs(e2); 

ypos=abs(z3(:,l)); 

ave=mean(elpos) 

aveft=mean(e2pos) 

outave=mean(ypos) 



% Average residual 
% Average residual (Prefiltered Data) 
% Average actuator output 
% Mean Square Fit calculated above 
% Mean Square Fit (Prefiltered) 



fit 

fitf 



% Calculate and plot actual output minus model output 

z3out=z3(:,l); 

for i=5000:5100; 

difft(i)=z3out(i)-yhf(i); 

end 

z3out=z3(:,l); 
for j=5000:5100; 
diff(j)=z3out(j)-yh(j); 



Figure(12) 

plot(diff,'r') 

axis([5000, 5 100,-0.5,0.5]) 
title('Actual Output - Model Output') 

Figure(13) 

plot(difft,'r') 

axis([5000, 5 100,-0.5,0.5]) 

title('Actual Output - Model Output *Data Filtered*') 
avediff=mean(abs(diff)) % Average output difference 

avedft=mean(abs(difft)) % Average output difference(Prefiltered) 

% Generate Zero-Pole Plot 
Figure( 14) 

zpplot 1 (th2zp(arxm),3 , [] , 1 .2) 

title([’Zero Pole Plot Box- Jenkins Model [',NB,' ',NC,' ',ND,' ',NF,' ',NK,'] ’]) 

Figure(15) 

zpplot 1 (th2zp(arxmft),3, [] ,1.2) 

title(['Zero Pole Plot Box-Jenkins Model [',NB,' ',NC,' ',ND,' ',NF,' ',NK,'] *Data 
Prefiltered*']) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% plotarx.m % 

% Compare Three ARX Model Structures % 

% George D. Beavers % 

% June 1997 % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



end 
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% This program generates the frequency response plots of the three best ARX model % 
% structures for evaluation and structure selection % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



j=input('enter strut»'); 

eval(['load 4h',int2str(j)]) 
eval(['load s',int2str(j),'m']) 

z2=[trace_y(j , 1 : 10000)' trace_y(7, 1 : 10000)'] ; 
z2=dtrend(z2); 

z3=[trace_y(j, 10001 :20000)' trace_y(7, 10001 :20000)']; 
z3=dtrend(z3); 

nnl=[7 6 1] % Top three structures 

nn2=[15 15 1] 
nn3=[12 12 1] 

Xl=o2ilx; 

Y 1 =20.*log 1 0(abs(o2i 1 )); 



zf=idfilt(z2, 10,. 52); 

% Construct the models % 
fml=arx(zf,nnl); 
fm 1 =sett(fm 1 ,0. 1 e-3) ; 
gfml=th2ff(fml); 

fm2=arx(zf,nn2); 
fm2=sett(fm2,0. le-3); 
gfm2=th2ff(fm2) ; 



fm3=arx(zf,nn3); 
fm3=sett(fm3,0. le-3); 
gfm3=th2ff(fm3); 

% Generate comparison plots 
Figure(l) 

semiiogx(X 1 , Y 1 ,'b'), 

title(['Strut #',int2str(j),' Magnitude']) 

xlabel('Frequency (Hz)') 

ylabel('Magnitude(dB)') 

hold on 
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semilogx(gfml(:,l)/(2*pi),20*Iogl0(gfml(:,2)),'r'); 
semilogx(gfm2(:,l)/(2*pi),20*logl0(gfm2(:,2)),'m'); 
semilogx(gfm3(:,l)/(2*pi),20*logl0(gfm3(:,2)),’w’);,grid; 
hold off 

Figure(2); 

eval(['load s',int2str(j),'p']) 

X2=o2ilx; 

Y2= 1 80*unwrap(angle(o2i 1 ))/pi; 
semilogx(X2,Y2,'b'), 
title(['Strut #',int2str(j),' Phase’]) 
xlabel(’Frequency (Hz)'); 
ylabel(’Deg'); 
hold on 

semilogx(gfm 1 (:, 1 )/(2*pi),(gfm 1 (:,4)),Y); 

semilogx(gfm2(:4)/(2*pi),(gfm2(:,4)),'m'); 

semilogx(gfm3(:,l)/(2*pi),(gfm3(:,4)),’w');,grid; 

hold off 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



% plotamb.m % 

% Model Comparison % 

% George D. Beavers % 

% June 1997 % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% This program generates the frequency response plots of the three best models based % 
% on the ARX structure^ 6 1]. ARX, ARMAX and BOX- JENKINS models are % 
% compared to select the best model % 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
j=input('enter strut»’); 

eval(['load 4h’,int2str(j)]) 
eval([’load s',int2str(j),'m’]) 

z2=[trace_y(j, 1:10000)' trace_y(7,l: 10000)']; 
z2=dtrend(z2); 

z3=[trace_y(j, 1 000 1 :20000)' trace_y(7, 10001 :20000)']; 
z3=dtrend(z3); 

% Define model structure % 
nnl=[7 6 2 1] 
nn2=[7 6 1] 
nn3=[6 2 2 7 1] 

Xl=o2ilx; 
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Y 1 =20. *log 1 0(abs(o2i 1 )); 

% Generate models (non-filtered) % 
%ml=armax(z2,nn 1 ,60); 
%ml=sett(ml,0.1e-3); 

%gm 1 =th2ff(m 1 ) ; 

%m2=arx(z2,nn2); 

%m2=sett(m2,0. le-3); 
%gm2=th2ff(m2) ; 
%m3=bj(z2,nn3,60); 

%m3=sett(m3,0. le-3); 
%gm3=th2ff(m3); 



zf=idfilt(z2 , 1 0,.52); 

% Generate Models (Prefiltered) 
fml=armax(zf,nnl ,60); 
fml=sett(fml ,0. le-3); 
gfml=th2ff(fml); 

fm2=arx(zf,nn2); 

fm2=sett(fm2,0.1e-3); 

gfm2=th2ff(fm2); 



fm3=bj(zf,nn3,60); 
fm3=sett(fm3,0. le-3); 
gfm3=th2ff(fm3); 

% Generate plots 
Figure(l) 

semilogx(X 1 , Y1 ,’b') 

title(['Strut #',int2str(j),' Magnitude']) 

xlabel (’Frequency (Hz)') 

y label( 'Magnitude(dB )') 

hold on 

semilogx(gfml(:,l)/(2*pi),20*logl0(gfml(:,2)),'r'); 
semilogx(gfm2(:,l)/(2*pi),20*logl0(gfm2(:,2)),’m'); 
semilogx(gfm3(:,l)/(2*pi),20*logl0(gfm3(:,2)),'w');,grid; 
hold off 

Figure(2); 

eval(['load s’,int2str(j),’p']) 

X2=o2ilx; 

Y2=180*unwrap(angle(o2il))/pi; 
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semilogx(X2,Y2,'b') 
title(['Strut #',int2str(j),' Phase']) 
xlabel('Frequency (Hz)'); 
ylabel('Deg'); 
hold on 

semilogx(gfm 1 (:, 1 )/(2*pi),(gfm 1 (:,4)),'r’); 
semilogx(gfm2(: , 1 )/(2*pi),(gfm2( : ,4)),'m'); 
semilogx(gfm3(: , 1 )/(2*pi),(gfm3(: ,4)),'w');,grid; 

hold off 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



% compid.m % 

% Interaction Comparison % 

% June 1997 % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% This program generates the frequency response plots of the actuators when actuator % 
% number one is excited with 50mV “pink noise”. These plots use the ARX[7 6 1] % 

% model frequency repines of the respective actuator. This is compared with the % 
% modeled response of actuator number one These plots are used to check for % 

°7o interaction between actuators % 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



function []=compid(j) 

% j= actuator number% 
int2str(j) 

disp(['load 4h' int2str(j)]) 
eval(['load ’ '4h' int2str(j) '.mat']) 

z 1 =[trace_y( 1,1:1 0000)' trace_y(7, 1 : 1 0000)’] ; 
zl=dtrend(zl); 

zl a=[trace_y( 1 , 10001 :20000)’ trace_y(7, 1 000 1 :20000)']; 
zla=dtrend(zla); 

thl=arx(zl,[7 6 1]); 
th 1 =sett(th 1,0.1 e-3) ; 
gth 1 =th2ff(th 1); 



for i=2:6 



z2=[trace_y(i, 1 : 1 0000)' trace_y(7, 1 : 1 0000)'] ; 
z2=dtrend(z2); 

z7=[trace_y(i, 1000 1:20000)' trace_y(7, 10001 :20000)']; 
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z7=dtrend(z7); 



th=arx(z2,[7 6 1]); 
th=sett(th,0.1e-3); 
gth=th2ff(th); 



Figure(i-l); 

semilogx(gth 1 (:, 1 )/(2*pi),20*log 10(gth 1 (:,2)),'b’); 
hold on 

semilogx(gth(:,l)/(2*pi),20*logl0(gth(:,2)),'r'); 
title('Magnitude: Input=50mV Pink Noise') 
xlabel('Frequency (Hz)') 
ylabel('dB') 

text( 10,45, ['Excite Actuator#' num2str(j) 'Blue; Sense Actuator ', 
num2str(i),'=Red']) 
grid 

hold off 



end 



for i=2:6 

z2=[trace_y (i, 1 : 1 0000)' trace_y(7, 1 : 1 0000)'] ; 
z2=dtrend(z2); 

z7=[trace_y(i, 10001 ;20000)' trace_y(7, 10001 :20000)’]; 
z7=dtrend(z7); 

th=arx(z2,[7 6 1]); 
th=sett(th,0.1e-3); 
gth=th2ff(th); 

Figure(i+4) 

semilogx(gthl(:,l)/(2*pi),(gthl(:,4)),'b'); 
hold on 

semilogx(gth(: , 1 )/(2*pi),(gth(:,4)),'r’); 
title('Phase: Input=50mv Pink Noise') 

text( 10,-550, ['Excite Actuator #' num2str(j) '=Blue; Sense Actuator ',num2str(i), 
'=Red']) 

xlabel('Frequency (Hz)') 

ylabel('Degrees') 

grid 

hold off 
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end 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% figl.m, fig2.m, fig3.m, fig4.m % 

% PSD of Applied Controllers % 

% June 1997 % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% This program generates the power spectral density plots resulting from applying % 
% a lead-lag compensator to an individual actuator and the all six actuators to evaluate % 
% the effects of actuator interaction % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%figl.m% 

% Plots over range of 0-5kHz% 
clear 

load s 1 %Open loop data with 42 Hz disturbance on% 

Figure(l) 

clg 

psd(trace_y(2, 1:20000), 1024*1, 10000, []); 

% [popen,fopen]=psd(trace_y(2, 1 :20000), 1 024* 1 2, 1 0000,[]); 
hold on 

load s2 %Compensator applied to actuator #2 ( 

psd(trace_y (2, 1 :20000), 1024* 1 , 1 0000,kaiser(5 1 2,5), 256); 

%[ps2,fs2]=psd(trace_y(2, 1:20000), 1024* 12, 10000, []); 

load s3 %Compensator applied to all six actuators 

psd(trace_y(2, 1 :20000),1024* 1 ,10000,[]); 

%[ps3,fs3]=psd(trace_y(2, 1:20000), 1024*12, 10000, []); 
hold off 

%semilogx((fopenk),20*logl0(popenk),'r',(fs2k),20*logl0(fs2k),y,(fs3k),20*logl0(fs3k) 

,'b') 

Title('Blue-Open Loop; Magenta-Strut #2 Closed Loop; Red-All 6 Struts Closed Loop' ) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%fig2.m% 

% Plots same Figure as figl.m over range of 0-200 Hz% 
clear 

load si %Open loop-42 Hz Disturbance on% 

Figure(2) 

clg 

psdb(trace_y(2, 1:20000), 1024*8, 10000, []); 

%[popen,fopen]=psd(trace_y(2, 1 :20000), 1 024* 1 2, 1 0000,[]); 
hold on 
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load s2 % Compensator applied to actuator #2% 

psd(trace_y(2,l:20000), 1024*8,10000, []); 

%[ps2,fs2]=psd(trace_y(2, 1:20000), 1024* 12,10000,[]); 

load s3 % Compensator applied to actuator #2% 

psd(trace_y(2, 1:20000), 1024*8, 10000,[]); 

%[ps3,fs3]=psd(trace_y(2, 1:20000), 1024* 12,10000,0); 
hold off 

%semilogx((fopenk),20*logl0(popenk),'r',(fs2k),20*logl0(fs2k),'y',(fs3k),20*logl0(fs3k) 

,'b') 

Title('Blue-Open Loop; Magenta-Strut #2 Closed Loop; Red-All 6 Struts Closed Loop' ) 
%Title('Open Loop-42Hz Disturbance' ) 
axis([0 200 -70 30]) 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%fig3.m% 

%plots PSD sensed by 3 geophones; actuator #2, #3 and #6 over range of 0-200Hz% 
clear 

load s2 % Compensator applied to actuator #2 % 

Figure(3) 

clg 

psd(trace_y(2, 1 :20000), 1 024* 1 , 1 0000, []); 

%[popen,fopen]=psd(trace_y(2,l :20000), 1024* 12, 10000, []); 
hold on 

psd(trace_y(3,l :20000), 1024* 1,10000, []); 

% [ps2,fs2]=psd(trace_y(3, 1 :20000), 1 024* 12,1 0000, []); 
psd(trace_y(6, 1 :20000), 1 024* 1 , 10000, []); 

%[ps3, fs3]=psd(trace_y(6, 1:20000), 1024* 12, 10000, []); 
hold off 

%semilogx((fopenk),20*logl0(popenk),'r',(fs2k),20*logl0(fs2k),'y',(fs3k),20*logl0(fs3k) 

;b’) 

Title('Blue-#2 Geophone; Magenta-#3 Geophone; Red-# 6 Geophone' ) 

%axis([0 200 -70 30]) 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%fig3.m% 

%plots PSD sensed by 3 geophones; actuator #2, #3 and #6 over range of 0-5kHz% 
clear 

load s2 %Compensator applied to actuator #2 % 

Figure(4) 

clg 

psd(trace_y(2, 1:20000), 1024*8, 10000, []); 

%[popen,fopen]=psd(trace_y(2, 1:20000), 1024* 1 2, 10000,[]); 
hold on 
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psd(trace_y(3,l :20000), 1024*8, 10000, []); 

%[ps2,fs2]=psd(trace_y(3,l :20000), 1024* 1 2, 10000, []); 
psd(trace_y(6,l :20000), 1024*8, 10000, []); 

% [ps3,fs3]=psd(trace_y(6, 1 :20000), 1 024* 1 2, 1 0000, []) ; 
hold off 

%semilogx((fopenk),20*logl0(popenk),'r',(fs2k),20*logl0(fs2k),y,(fs3k),20*logl0(fs3k) 

,V) 

Title('Blue-#2 Geophone; Magenta-#3 Geophone; Red-# 6 Geophone' ) 
axis([0 200 -70 30]) 



ill 



APPENDIX B. MODEL STRUCUTURE SELECTION PLOTS 



A. Structure Acceptance Plots for Actuator Two 

The following plots were utilized in selecting the ARX structure [na=7 nb=6 nk=l] 
for actuator number two. 

RETURN TO COMMAND SCREEN TO SELECT # OF PARAMETERS TO BE ESTIMATED 
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Figure A. 1 

Parameter Number Vs. Loss Function For nk= 1 
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Zero Pole Plot ARX Model [7 6 1] Model ‘Data Prefiltered* 




Figure A. 2 
Zero-Pole Plot 
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Phase (Deg) , Magnitude(dB) 




Frequency (Hz) 



(a) Magnitude 




Frequency (Hz) 

(b) Phase 
Figure A.3 

Frequency Response Comparison 
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Output Magnitude 



Plot of ARX MODEL [7 6 1] Output vs Actual output *Data Prefiltered* 




Sample 
Figure A.4 

Model Output Vs. Actual Output 



Correlation function of residuals. Output # 1 




Cross corr. function between input 1 and residuals from output 1 




Figure A. 5 

Auto and Cross Correlation Functions 
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Residual 




5000 5020 5040 5060 5080 5100 

Sample 

Figure A. 6 
Residual Vs. Sample 



Actual Output - ARX Model [7 6 1] Model *Data Prefiltered* 




Figure A. 7 

Actual Output and Model Output Difference 
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B. Model Structure Acceptance for Actuator Three 

The following plots were utilized in selecting the ARX structure [na=7 nb=7 nk=l] 
for actuator number three. 



RETURN TO COMMAND SCREEN TO SELECT # OF PARAMETERS TO BE ESTIMATED 
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Zero Pole Plot ARX Model [7 6 1] Model 




Figure B.2 
Zero-Pole Plot 
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Frequency (Hz) 
(a) Magnitude 




Frequency (Hz) 
(b) Phase 



Figure B.3 

Frequency Response Comparison 
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Output Magnitude 



Plot of ARX MODEL [761] Output vs. Actual output 




Sample 
Figure B.4 

Model Output Vs Actual Output 



Correlation function of residuals. Output# 1 




Cross corr. function between input 1 and residuals from output 1 




Figure B.5 

Auto and Cross Correlation Functions 
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Residual of ARX Model [7 6 1] 




Figure B .6 
Residual Vs Sample 



Actual Output - ARX Model [7 6 1] Model 




Figure B.7 

Actual Output and Model Output Difference 
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C. Model Structure Acceptance for Actuator Four 

The following plots were utilized in selecting the ARX structure [na=7 nb=7 nk=l] 
for actuator number four. 



RETURN TO COMMAND SCREEN TO SELECT ff OF PARAMETERS TO BE ESTIMATED 




Figure C. 1 

Parameter Number Vs. Loss Function For nk= 1 



Zero Pole Plot ARX Model [7 61] Model ‘Data Prefiltered* 




Figure C.2 
Zero-Pole Plot 
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Frequency (Hz) 



(a) Magnitude 




Frequency (Hz) 
(b) Phase 



Figure C.3 

Frequency Response Comparison 
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Output Magnitude 



Plot of ARX MODEL [7 6 1] Output vs. Actual output ‘Data Prefiltered* 




Sample 



Figure C.4 

Model output Vs. Actual Output 



Correlation function of residuals. Output# 1 




Cross conr. function between input 1 and residuals from output 1 




Figure C.5 

Auto and Cross Correlation Functions 



123 



Residual 




5000 5020 5040 5060 5080 5100 

Sample 

Figure C.6 
Residual Vs. Sample 

Actual Output - ARX Model [7 6 1] Model ‘Data Prefiltered* 




Figure C.7 

Actual Output and Model Output Difference 
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D. Model Structure Acceptance for Actuator Five 

The following plots were utilized in selecting the ARX structure [na=7 nb=7 nk=T] 
for actuator number five. 

RETURN TO COMMAND SCREEN TO SELECT # OF PARAMETERS TO BE ESTIMATED 




Figure D.l 

Parameter Number Vs. Loss Function For nh= 1) 



Zero Pole Plot ARX Model [7 6 1] Model 'Data Prefiltered* 
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Frequency (Hz) 
(a) Magnitude 




Frequency (Hz) 

(b) Phase 
Figure D.3 

Frequency Response Comparison 
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Output Magnitude 



Plot of ARX MODEL [7 6 1] Output vs. Actual output ‘Data Prefiltered* 




Sample 



Figure D.4 

Model Output Vs Actual Output 



Correlation function of residuals. Output # 1 




Cross corr. function between input 1 and residuals from output 1 




Figure D.5 

Auto and Cross Correlation Functions 
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Residual 




5000 5020 5040 5060 5080 5100 

Sample 

Figure D.6 
Residual Vs. Sample 



Actual Output - ARX Model [7 6 1] Model *Data Prefiltered* 




Figure D.7 

Actual Output and Model Output Difference 
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E. Model Structure Acceptance for Actuator Six 

The following plots were utilized in selecting the ARX structure [na=7 nb=7 nk=l] 
for actuator number six. 

RETURN TO COMMAND SCREEN TO SELECT # OF PARAMETERS TO BE ESTIMATED 




Figure E. 1 

Parameter Number Vs. Loss Function For nk= 1 




Figure E.2 
Zero-Pole Plot 
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Frequency (Hz) 
(a) Magnitude 




Frequency (Hz) 

(b) Phase 
Figure E.3 

Frequency Response Comparison 
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Output Magnitude 



Plot of ARX MODEL [7 6 1] Output vs. Actual output ’Data Prefiltered* 




Figure E.4 

Model Output Vs. Actual Output 



Correlation function of residuals. Output # 1 




Cross corr. function between input 1 and residuals from output 1 




Figure E.5 

Auto and Cross Correlation Functions 
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Residual 
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Figure E.6 
Residual Vs Sample 

Actual Output - ARX Model [7 6 1] Model *Data Prefiltered* 




Figure E.7 

Actual Output and Model Output Difference 
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APPENDIX C. MODEL VALIDATION PLOTS 



A. Application of ARMAX Method to Actuator Number One 

The following plots were used in evaluating the application of the ARMAX and 
Box-Jenkins parameter estimation methods in modeling actuator number one. The 
structure [na=7 nb=6 nc=2 nk=l] was used. 



Blue=Actuator#1 Red=ARMAX MODEL[7 6 21] “Data Prefiltered* 




Frequency (Hz) 

(a) Magnitude 

Blue=Actuator#1 Red=ARMAX MODEL[7 6 2 1] ‘Data Prefiltered* 




Frequency (Hz) 



(b) Phase 
Figure A. 1 
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Frequency Response Comparison 

Plot of ARMAX MODEL[7 6 2 1] Output vs. Actual output ‘Data Prefiltered* 




Sample 



Figure A.2 

Model Output Vs. Actual Output 



Correlation function of residuals. Output # 1 




Cross corn function between input 1 and residuals from output 1 




Figure A. 3. 

Auto and Cross-Correlation Functions 
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Residual 
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Figure A.4 
Residual Vs. Sample 



Actual Output- ARMAX[7 6 2 IJModel Output ‘Data Filtered* 




Figure A. 5 

Actual Output and Model Output Difference 
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Zero Pole Plot ARMAX [7 6 2 1] Model ‘Data Prefiltered* 
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Figure A 6 
Zero-Pole Plot 
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B. Application of ARMAX Method to Actuator Number Two 



The following plots were used in evaluating the application of the ARMAX 
parameter estimation method in modeling actuator number two. The structure [na=7 nb=6 
nc=2 nk=l] was used. Two additional plots from the Box-Jenkins \nb=6 nc = 2 /?c= 2 nci=l 
nk= 1] model are included for comparison. 



Blue=Actuator#2Red=ARMAX MODEL[7 6 2 1] *Data Prefiltered* 




Frequency (Hz) 



(a) Magnitude 



Blue=Actuator#2Red=ARMAX MODEl[7 6 2 1] ‘Data Prefiltered* 




Frequency (Hz) 



(b) Phase 



137 



Figure B.l 

Frequency Response Comparison 
Plot of ARMAX MODEL[7 6 2 1] Output vs. Actual output 'Data Prefiltered' 




Figure B.2 

Model Output Vs. Actual Output 



Correlation function of residuals. Output# 1 




Cross con. function between input 1 and residuals from output 1 




Figure B.3 

Auto and Cross Correlation Functions 
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Residual 
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Figure B.4 
Residual Vs Sample 

Actual Output - ARMAX[7 6 2 1]Model Output *Data Filtered* 




Figure B.5 

Actual Output and Model Output Difference 
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Zero Pole Plot ARMAX [7 6 2 1] Model *Data Prefiltered 
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Figure B.6 
Zero-Pole Plot 

Zero Pole Plot Box-Jenkins Model [6227 1] ‘Data Prefiltered* 




Figure B.7 

Zero-Pole Plot for Box-Jenkins Model 
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Correlation function of residuals. Output# 1 




Cross corr. function between input 1 and residuals from output 1 




lag 

Figure B.8 

Auto and Cross Correlation Functions for Box-Jenkins Model 
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C. Application of ARMAX Method to Actuator Number Three 



The following plots were used in evaluating the application of the ARMAX 
parameter estimation method in modeling actuator number three. The ARX based 
structure [na=7 nb=6 nc=2 nk=l] was used. 



Blue=Actuator#3Red=ARMAX M0DEL[7 6 21] ‘Data Prefiltered* 




Frequency (Hz) 



(a) Magnitude 

Blue=Actuator#3Red=ARMAX MODEL[7 6 21] *Data Prefiltered* 




(b) Phase 
Figure C.l 

Frequency Response Comparison 
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Plot of ARMAX MODEL[7 6 2 1] Output vs. Actual output ‘Data Prefiltered* 




Sample 



Figure C.2 

Model Output Vs. Actual Output 



Correlation function of residuals. Output# 1 




Cross corr. function between input 1 and residuals from output 1 




Figure C.3 

Auto and Cross Correlation Functions 
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Residua! 
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Figure C.4 
Residual Vs. Sample 

Actual Output - ARMAX[7 6 2 1 ]Model Output *Data Filtered* 




Figure C.5 

Actual Output and Model Output Difference 
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Zero Pole Plot ARMAX [7 6 2 1] Model *Data Prefiltered* 
1.5 

1 

0.5 
0 

-0.5 
-1 
-1.5 

-1.5 -1 -0.5 0 0.5 1 1.5 

Figure C.6 
Zero-Pole Plot 
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D. Application of ARMAX Method to Actuator Number Four. 



The following plots were used in evaluating the application of the ARMAX 
parameter estimation method in modeling actuator number four. The structure [na=7 nb=6 
nc=2 nk=l] was used. 



Blue=Actuator#4Red=ARMAX MODEL[7 6 2 1] 

60 1 i-- i | mm | | 1 1 m i r 




10 10 10 10 10 
Frequency (Hz) 

(a) Magnitude 




Frequency (Hz) 
(b) Phase 
Figure D 1 
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Frequency Response Comparison 
Plot of ARMAX MODEL[7 6 2 1] Output vs. Actual output 




Figure D.2 

Model Output Vs. Actual Output 



Correlation function of residuals. Output # 1 




Cross corr. function between input 1 and residuals from output 1 




Figure D.3 

Auto and Cross Correlation Functions 
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Residual 
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Figure D.4 
Residual Vs. Sample 



Actual Output - ARMAX[7 6 21] Model Output 




Figure D.5 

Actual Output and Model Output Difference 
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Zero Pole Plot ARMAX M0DEL[7 6 2 1] 
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Figure D.6 
Zero-Pole Plot 



149 



E. Application of ARMAX Method to Actuator Number Five. 



The following plots were used in evaluating the application of the ARMAX 
parameter estimation method in modeling actuator number five. The structure [na=7 nb=6 
nc=2 nk=l] was used. 




Frequency (Hz) 



(a) Magnitude 




Frequency (Hz) 



(b) Phase 
Figure E. 1 

Frequency Response Comparison 
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Plot of ARMAX MODEL[7 6 2 1] Output vs. Actual output ‘Data Pre filtered* 




Sample 



Figure E.2 

Model Output Vs. Actual Output 



Correlation function of residuals. Output # 1 




Cross corn function between input 1 and residuals from output 1 




lag 

Figure E.3 

Auto and Cross Correlation Functions 



151 



Residual 
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Figure E.4 
Residual Vs. Sample 



Actual Output - ARMAX[7 6 2 1]Model Output *Data Filtered* 




Figure E.5 

Actual Output and Model Output Difference 
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Zero Pole Plot ARMAX [7 6 2 1] Model ‘Data Prefiltered* 
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Figure E.6 
Zero-Pole Plot 
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F. Application of ARJVIAX Method to Actuator Number Six. 

The following plots were used in evaluating the application of the ARMAX 
parameter estimation method in modeling actuator number six. The structure[na=7 nb=6 
nc=2 nk=l] was used. 



Blue=Actuator#6Red=ARMAX MODEL[7 6 2 1] *Data Prefiltered* 
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Blue=Actuator#6Red=ARMAX MODEL[7 6 2 1] *Data Prefiltered* 
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(b) Phase 
Figure F. 1 
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Frequency Response Comparison 

Plot ofARMAX MODEL[7 6 2 1] Output vs. Actual output ‘Data Prefiltered* 




Figure F.2 

Model Output Vs. Actual Output 



Correlation function of residuals. Output # 1 




Cross com function between input 1 and residuals from output 1 




Figure F.3 

Auto and Cross Correlation Functions 
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Residual 



Residual of ARMAX MODEL[7 6 2 1] ‘Data Prefiltered* 
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Figure F.5 

Actual Output and Model Output Difference 
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Zero Pole Plot ARMAX [762 1] Model *Data Prefiltered* 
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Figure F.6 
Zero Pole Plot 
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